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Abstract 

We  consider  an  electromagnetic  interrogation  technique  in  two  and  three  dimensions  for  identi¬ 
fying  the  dielectric  parameters  (including  the  permittivity,  the  conductivity  and  the  relaxation  time) 
of  a  Debye  medium.  In  this  technique  a  travelling  acoustic  pressure  wave  in  the  Debye  medium  is 
used  as  a  virtual  reflector  for  an  interrogating  microwave  electromagnetic  pulse  that  is  generated 
in  free  space.  The  reflections  of  the  microwave  pulse  from  the  air-Debye  interface  and  from  the 
acoustic  pressure  wave  are  recorded  at  a  remote  antenna.  The  data  is  used  in  an  inverse  problem  to 
estimate  the  locally  pressure  dependent  dielectric  parameters  of  the  Debye  medium.  We  present  a 
time  domain  formulation  that  is  solved  using  finite  differences  (FDTD)  in  time  and  in  space.  Per¬ 
fectly  matched  layer  (PML)  absorbing  boundary  conditions  are  used  to  absorb  outgoing  waves  at 
the  finite  boundaries  of  the  computational  domain,  preventing  spurious  reflections  from  reentering 
the  domain. 

Using  the  method  of  least  squares  for  the  parameter  identification  problem,  we  compare  two  dif¬ 
ferent  algorithms  (the  gradient  based  Levenberg-Marquardt  method,  and  the  gradient  free,  simplex 
based  Nelder-Mead  method)  in  solving  an  inverse  problem  to  calculate  estimates  for  two  or  more 
dielectric  parameters.  Finally  we  use  statistical  error  analysis  to  construct  confidence  intervals  for 
all  the  presented  estimates,  thereby  providing  a  probabilistic  statement  about  the  computational 
procedure  with  uncertainty  aspects  of  estimates. 

Keywords:  Electromagnetic-acoustic  interaction,  Debye  dielectric  materials,  pulsed  antenna 
source  microwaves,  inverse  problems,  FDTD,  statistical  inference. 
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1  Introduction 


The  study  of  transient  electromagetic  waves  in  lossy  dispersive  dielectrics  is  of  great  interest 
due  to  numerous  applications  [AMP94,  APM89].  These  include  microwave  imaging  for 
medical  applications  in  which  one  seeks  to  investigate  the  internal  structure  of  an  object 
(the  human  body)  by  means  of  electromagnetic  fields  at  microwave  frequencies  (300  MHz 
-  30  GHz).  This  is  done  for  example,  to  detect  cancer  or  other  anomalies  by  studying 
changes  in  the  dielectric  properties  (such  as  permittivity,  conductivity,  relaxation  times, 
etc.)  of  tissues  [FMS03].  One  generates  the  electrical  property  distributions  in  the  body 
(. Microwave  maps)  with  the  hope  that  such  properties  of  different  bodily  tissues  are  related 
to  their  physiological  state.  With  such  noninvasive  interrogating  techniques  one  can  study 
properties  and  defects  in  biological  tissues  with  very  little  discomfort  to  the  subjects.  Other 
potential  applications  for  such  interrogation  techniques  are  nondestructive  damage  detection 
in  aircraft  and  spacecraft  where  very  high  frequency  electromagnetic  pulses  can  be  used 
to  detect  the  location  and  width  of  cracks  that  may  be  present  [BGW03].  Additional 
applications  are  found  in  mine,  ordinance  and  camouflage  surveillance,  and  subsurface  and 
atmospheric  environmental  modelling  [BBL00]. 

Another  important  technique  employed  in  noninvasive  interrogation  of  objects  is  the  use 
of  ultrasonic  waves  in  industrial  and  medical  applications.  The  interaction  of  electromag¬ 
netic  and  sound  waves  in  matter  has  been  widely  studied  in  acoustooptics  [And67,  Kor97]. 
Brillouin  first  showed  that  electromagnetic  and  sound  waves  can  interact  in  a  medium 
and  influence  each  others  propagation  dynamics  [Bri60] .  In  particular,  the  response  of  the 
atomic  electrons  in  a  material  medium  to  the  applied  electric  field  results  in  a  polarization 
of  the  material  with  a  resulting  index  of  refraction  that  is  a  function  of  the  density  in  the 
material.  Consequently,  the  material  density  fluctuations  produced  by  a  sound  wave  induce 
perturbations  in  the  index  of  refraction.  Thus  an  electromagnetic  wave  transmitted  through 
the  medium  will  be  modulated  by  the  sound  wave,  and  scattering  and  reflection  will  occur. 
Conversely,  a  spatial  variation  of  the  electromagnetic  field  intensity  and  the  corresponding 
electromagnetic  stress  lead  to  a  volume  force  distribution  in  the  medium.  This  is  called 
electrostriction,  and  can  lead  to  sound  generation.  This  mutual  interaction  between  light 
and  sound  may  also  lead  to  instabilities  and  wave  amplification  [MI68]. 

The  work  in  this  paper  is  based  on  ideas  underlying  the  techniques  formulated  in  [BBL00] 
and  [ABR02],  In  [BBL00],  the  authors  present  an  electromagnetic  interrogation  technique 
for  general  dispersive  media  backed  by  either  a  perfect  conducting  wall  to  reflect  the  elec¬ 
tromagnetic  waves,  or  by  using  standing  acoustic  waves  as  virtual  reflectors.  In  [ABR02], 
this  work  was  extended  to  using  travelling  waves  as  acoustic  reflectors.  In  both  cases  the 
techniques  were  demonstrated  in  one  dimension  by  assuming  plane  electromagnetic  waves 
impinging  at  normal  incidence  on  slabs  of  dispersive  media.  A  one  dimensional  setting  with 
normally  incident  waves  is  not  geometrically  feasible  in  many  applications.  In  [BB04a]  a 
specialized  multidimensional  version  of  the  electromagnetic  interrogation  problem  is  pre¬ 
sented  for  dispersive  media  backed  by  perfect  conductor  walls  and  the  finite  difference  time 
domain  (FDTD)  method  is  utilized  [Taf95,  Taf98].  However,  the  object  of  interrogation  may 
not  be  backed  by  perfect  conductor  walls  in  many  important  applications,  and  hence  one 
seeks  another  type  of  interface  which  will  reflect  the  propagating  electromagnetic  pulses.  It 
is  therefore  of  great  importance  to  ascertain  (computationally  and  experimentally)  whether 
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Figure  1:  Antenna  and  dielectric  slab  configuration. 


the  interrogation  ideas  of  [ABR02]  using  acoustic  pressure  waves  as  virtual  reflectors  can 
be  effective  in  multidimensional  settings.  The  goal  of  this  paper  is  the  investigation  of  com¬ 
putational  and  statistical  aspects  of  this  methodology  in  higher  dimensions.  At  present  we 
are  also  trying  to  validate  these  proposed  techniques  by  testing  the  ideas  in  this  paper  in 
an  experimental  setup.  Our  group  has  built  an  antenna  that  we  are  using  to  study  the 
acoustooptical  interaction  in  agar-agar  in  which  the  sound  waves  are  produced  by  means 
of  a  transducer  [ABK04],  In  such  a  setup,  as  well  as  in  many  applications,  one  cannot  as¬ 
sume  normally  incident  plane  waves  impinging  on  the  agar,  and  hence  one  needs  to  consider 
the  more  general  case  of  oblique  incidence.  This  further  motivates  the  study  in  a  multi¬ 
dimensional  setting  of  an  electromagnetic  interrogation  technique  using  the  acoustooptical 
interaction  to  investigate  dielectic  properties  of  dispersive  media. 

We  begin  with  Maxwell’s  equations  for  a  first  order  dispersive  (Debye)  media  with 
Ohmic  conductivity  and  incorporate  a  model  that  describes  the  electromagnetic/acoustic 
interaction.  Our  model  features  modulation  of  the  material  polarization  by  the  pressure 
wave  and  thus  the  behavior  of  the  electromagnetic  pulse.  This  approach  is  based  on  ideas 
found  in  [MI68].  We  incorporate  a  pressure  dependent  Debye  model  for  orientational  po¬ 
larization  into  Maxwell’s  equations.  We  model  the  propagation  of  a  non-harmonic  pulse 
from  a  finite  antenna  source  in  free  space  impinging  obliquely  on  a  planar  interface  into 
the  Debye  medium.  We  use  the  finite  difference  time  domain  method  (FDTD)  [Taf95] 
to  discretize  Maxwell’s  equations  and  to  compute  the  components  of  the  electric  field  in 
the  case  where  the  signal  and  the  dielectric  parameters  are  independent  of  the  y  variable. 
Figure  1  depicts  the  Debye  dielectric  slab  geometry  and  the  antenna  that  we  model  in 
our  problem  with  the  infinite  dielectric  slab  perpendicular  to  the  2-axis,  uniform  in  the 
(/-direction  lying  in  zi  <  2  <  1,  and  the  antenna  finite  in  the  x-direction  lying  in  the  region 
—00  <  y  <  00,  x\  <  x  <  X2,  respectively.  An  alternating  current  along  the  x-direction 
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of  the  antenna  located  at  z  =  zc  then  produces  an  electric  field  that  is  uniform  in  y  with 
nontrivial  components  Ex  and  E~  depending  on  (t,  x,  z).  When  propagated  in  the  xz- plane, 
this  results  in  oblique  incident  waves  on  the  dielectric  surface  in  the  xy-plane  at  z  =  z\.  (In 
[BF95]  the  authors  use  Fourier-series  in  the  frequency  domain  to  compute  the  propagation 
of  a  time  harmonic  pulse  train  of  plane  waves  that  enter  a  dielectric  across  a  planar  bound¬ 
ary  at  oblique  incidence.  They  showed  that  precursors  are  excited  by  short  rise  time  pulses 
at  oblique  incidence.  However,  the  use  of  Fourier  series  restricts  this  approach  to  harmonic 
pulses.)  The  oblique  waves  impinge  on  the  planar  interface  that  separates  air  and  the  Debye 
medium  at  z  =  z\  where  they  undergo  partial  transmission  and  partial  reflection.  A  re¬ 
flected  wave  travels  back  to  the  antenna  at  z  =  zc  and  is  recorded.  The  transmitted  part  of 
the  wave  travels  through  the  Debye  medium  and  intercepts  an  acoustic  pressure  wave  that 
is  generated  in  the  dielectric  medium.  Here  the  microwave  pulse  undergoes  partial  transmis¬ 
sion  and  reflection.  The  reflected  wave  travels  back  through  the  Debye  medium,  undergoes 
a  secondary  reflection/transmission  at  the  air-Debye  interface  and  is  then  recorded  at  the 
antenna.  Thus  the  antenna  records  three  pieces  of  information;  the  input  source,  the  reflec¬ 
tion  of  the  electromagnetic  pulse  from  the  air-Debye  interface  and  the  reflection  from  the 
pressure  wave  region.  This  data  is  recorded  at  the  center  of  the  antenna  and  will  be  used 
in  an  inverse  problem  to  estimate  the  dielectric  properties  of  the  Debye  medium.  We  note 
that  the  acoustic  wave  speed  is  many  orders  of  magnitude  smaller  than  the  speed  of  light. 
Thus  the  acoustic  pressure  wave  is  almost  stationary  relative  to  the  microwave  pulse.  In 
this  paper,  we  have  assumed  that  the  acoustic  pressure  wave  is  specified  a  priori  and  hence 
we  assume  that  the  pressure  wave  is  not  modified  by  the  electromagnetic  wave.  This  leads 
to  a  much  simplified  electromagnetic/acoustic  interaction  model.  The  dynamic  interaction 
between  the  two  and  the  appropriateness  of  this  simplifying  assumption  are  currently  being 
investigated  and  will  be  the  topic  of  a  future  paper. 

An  important  computational  issue  in  modelling  two  or  three  dimensional  wave  propaga¬ 
tion  is  the  construction  of  a  finite  computational  domain  and  appropriate  boundary  condi¬ 
tions  for  terminating  this  domain  in  order  to  numerically  simulate  an  essentially  unbounded 
domain  problem.  We  surround  the  domain  of  interest  by  perfectly  matched  absorbing  lay¬ 
ers  (PML’s),  thus  producing  a  finite  computational  domain.  The  original  PML  formulation 
(called  the  split  field  PML)  is  due  to  J.  P.  Berenger  [Ber94]  and  was  based  on  a  split  form 
of  Maxwell’s  equations  in  which  all  the  fields  themselves  were  split  into  orthogonal  subcom¬ 
ponents.  Many  different  variations  of  the  PML  are  now  available  in  the  literature.  Most 
of  these  different  formulations,  though  equivalent  to  the  split  field  PML,  do  not  involve  the 
splitting  of  the  fields  or  splitting  of  Maxwell’s  equations.  The  PML  formulation  that  we  use 
in  this  paper  is  based  on  the  anisotropic  (uniaxial)  formulation  that  was  first  proposed  by 
Sacks  et.  al.  [SKLL95]. 

We  discuss  numerical  results  for  the  forward  problem  for  two  different  test  cases.  The 
relaxation  times  of  the  Debye  media  in  these  two  test  cases  differ  by  many  orders  of  magni¬ 
tude.  We  first  perform  a  sensitivity  analysis  in  each  test  case  to  determine  which  parameters 
are  most  likely  to  be  accurately  estimated  in  an  inverse  problem.  We  also  determine  which 
of  the  pressure  dependent  parameters  are  the  most  significant  in  each  of  the  test  cases.  As 
we  shall  see,  the  significant  difference  in  the  orders  of  magnitude  of  the  relaxation  times  for 
the  two  test  cases  produces  inverse  problems  with  very  different  estimation  properties. 

We  next  discuss  a  parameter  estimation  problem  in  which  we  estimate  model  parameters 
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for  the  two  different  Debye  media  from  simulated  data.  We  formulate  our  inverse  problem 
using  the  method  of  least  squares.  We  compare  two  different  algorithms  in  solving  the 
parameter  estimation  problems;  the  Levenberg-Marquardt  method,  which  is  a  gradient 
based  algorithm,  and  the  Nelder-Mead  method,  which  is  based  on  simplices  and  is  gradient 
free.  For  the  first  test  case,  these  two  algorithms  produce  similar  results.  However,  in  the 
second  test  case,  we  see  the  presence  of  many  local  minima  and  the  Nelder-Mead  algorithm 
converges  to  a  local  minima  instead  of  the  global  minimum. 

We  also  consider  uncertainty  in  estimates  due  to  data  uncertainty.  We  present  a  model 
for  generating  noisy  data  that  we  use  in  our  inverse  problems  to  simulate  experimental 
settings  in  which  noise  enters  into  the  problem  in  a  natural  fashion.  We  then  discuss 
a  statistical  error  methodology  to  analyze  the  results  of  the  inverse  problems  based  on 
noisy  data.  This  statistical  error  analysis  is  used  to  obtain  confidence  intervals  for  all  the 
estimated  parameters,  thereby  providing  levels  of  confidence  that  can  be  associated  with 
estimates  obtained  with  our  methodology. 


2  Model  formulation:  The  forward  problem 


We  consider  Maxwell’s  equations  which  govern  the  electric  field  E  and  the  magnetic  field 
H  in  a  domain  H  with  charge  density  p.  Thus  we  consider  the  system 


<9D 

(i)  —  +J-  VxH  =  0,  in  (0,T)  x  fi, 
<9B 

(ii)  -7^-  +  V  x  E  =  0,  in  (0,  T)  x  P, 

(iii)  V  •  D  =  p,  in  (0,  T)  x  H, 

(iv)  V  •  B  =  0,  in  (0,  T)  x  H, 

(v)  E  x  n  =  0,  on  (0,  T)  x  dQ, 

(vi)  E(0,x)  =  0,  H(0,x)  =  0,  in  H. 


(1) 


The  fields  D,  B  are  the  electric  and  magnetic  flux  densities  respectively.  All  the  fields  in 
(1)  are  functions  of  position  r  =  (x,  y,  z)  and  time  t.  We  have  J  =  Jc  +  Js,  where  Jc  is  a 
conduction  current  density  and  Js  is  the  source  current  density.  We  assume  only  free  space 
(actually  the  antenna  in  our  example)  can  have  a  source  current,  and  Jc  is  only  found  in 
the  dielectric  material.  The  condition  (1,  v),  where  n  is  the  unit  outward  normal  to  d H, 
is  a  perfectly  conducting  boundary  condition  on  the  computational  domain  which  we  will 
replace  with  absorbing  conditions  in  the  sequel.  With  perfect  conductor  conditions,  the 
computational  boundary  acts  like  a  hard  wall  to  impinging  electromagnetic  waves  causing 
spurious  reflections  that  would  contaminate  the  data  that  we  wish  to  use  in  our  inverse 
problems. 

Constitutive  relations  which  relate  the  electric  and  magnetic  fluxes  D ,  B  and  the  electric 
current  Jc  to  the  electric  and  magnetic  fields  E,  H  are  added  to  these  equations  to  make  the 
system  fully  determined  and  to  describe  the  response  of  a  material  to  the  electromagnetic 
fields.  In  free  space,  these  constitutive  relations  are  D  =  eoE,  and  B  =  /ioH,  where  eo  and 
Ho  are  the  permittivity  and  the  permeability  of  free  space,  respectively,  and  are  constant 
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[Jac99] .  In  general  there  are  different  possible  forms  for  these  constitutive  relationships.  In  a 
frequency  domain  formulation  of  Maxwell’s  equations,  these  are  usually  converted  to  linear 
relationships  between  the  dependent  and  independent  quantities  with  frequency  dependent 
coefficient  parameters.  We  will  consider  the  case  of  a  dielectric  in  which  magnetic  effects  are 
negligible,  and  we  assume  that  Ohm’s  law  governs  the  electric  conductivity.  Thus,  within 
the  dielectric  medium  we  have  constitutive  relations  that  relate  the  flux  densities  D,B  to 
the  electric  and  magnetic  fields,  respectively.  We  have 

'  (i)  D  =e0E  +  P ID, 

(ii)  B  =moH,  (2) 

,  (hi)  Jc  =  ctE IDi 

In  (2),  Id  denotes  the  indicator  function  on  the  Debye  medium.  Thus,  Jc  =  0  in  air.  The 
quantity  P  is  called  the  electric  polarization.  It’s  contribution  to  the  electric  flux  is  zero 
in  air  and  is  nonzero  in  the  dielectric.  Electric  polarization  may  be  defined  as  the  electric 
field  induced  disturbance  of  the  charge  distribution  in  a  region.  This  polarization  may  have 
an  instantaneous  component  as  well  as  ones  that  do  not  occur  instantaneously;  the  latter 
usually  have  associated  time  constants  called  the  relaxation  times  which  are  denoted  by  r. 
We  let  the  instantaneous  component  of  the  polarization  to  be  related  to  the  electric  field 
by  means  of  a  dielectric  constant  eoX  and  denote  the  remainder  of  the  electric  polarization 
as  Pr.  We  have 

P  =  Pi  +  Pr  =  eoxE  +  Pr,  (3) 

and  hence  the  constitutive  law  (2,  i)  becomes 

D  =  eoer¥j  +  Pr,  (4) 


where  er  =  (1  +  x)  is  the  relative  permittivity  of  the  dielectric  medium. 

To  describe  the  behavior  of  the  media’s  macroscopic  electric  polarization  Pr,  we  employ 
a  general  integral  equation  model  in  which  the  polarization  explicitly  depends  on  the  past 
history  of  the  electric  field.  This  model  is  sufficiently  general  to  include  microscopic  po¬ 
larization  mechanisms  such  as  dipole  or  orientational  polarization  (Debye),  as  well  as  ionic 
and  electronic  polarization  (Lorentz),  and  other  frequency  dependent  polarization  mecha¬ 
nisms  [And67]  as  well  as  multiples  of  such  mechanisms  represented  by  a  distribution  of  the 
associated  time  constants  (e.g.,  see  [BG04]).  The  resulting  constitutive  law  can  be  given  in 
terms  of  a  polarization  or  displacement  susceptibility  kernel  g  as 


Pr  (t,x,z) 


s,  x,  z)E(s,  x,  z)ds, 


(5) 


which,  in  the  case  (Debye)  of  interest  to  us  here,  can  be  written  as 


g(t)  =e~^e°(e^  £oo). 

Such  a  Debye  model  can  be  represented  in  differential  form  as 


(6) 


tPr  +  Pr  =  e0(es  -  e^E, 
D  =  eoCooE  +  Pr, 


(7) 
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inside  the  dielectric,  whereas  Pr  =  0,  =  1  in  air,  and  =  er  inside  the  dielectric.  We 

will  henceforth  denote  Pr  by  P.  In  equations  (5)  and  (7),  the  parameters  es  and  denote 
the  static  relative  permittivity,  and  the  value  of  permittivity  for  an  extremely  high(~  oo) 
frequency  field,  respectively.  The  constant  £oo  is  called  the  infinite  frequency  permittiv¬ 
ity.  Biological  cells  and  tissues  display  very  high  values  of  the  dielectric  constants  at  low 
frequencies,  and  these  values  decrease  in  almost  distinct  steps  as  the  excitation  frequency 
is  increased.  This  frequency  dependence  of  permittivity  is  called  dielectric  dispersion  and 
permits  identification  and  investigation  of  a  number  of  completely  different  underlying  mech¬ 
anisms.  This  property  of  dielectric  materials  makes  the  problem  of  parameter  identification 
an  important  as  well  as  very  useful  topic  for  investigation. 

We  will  modify  system  (1)  and  the  constitutive  laws  (2)  by  performing  a  change  of 
variables  that  renders  the  system  in  a  form  that  is  convenient  for  analysis  and  computation. 
From  (1,  i)  we  have 

—  +  J  Jc(s,x)  ds^  —  VxH  =  —  Js,  in  (0 ,  T)  x  Q.  (8) 


Next,  we  define  the  new  variable 

D(f,  x)  =  D(f,  x)  +  f  Jc(s,x)ds. 

Jo 

Using  definition  (9)  in  (8)  and  (1)  we  obtain  the  modified  system 

<9D 

(i)  —  -VxH  =  -J,in(0,T)xU, 
<9B 

(ii)  +  VxE  =  0  in  (0, T )  x  U, 

(iii)  V  ■  D  =  0  in  (0,  T)  x  fi, 

(iv)  V  •  B  =  0  in  (0,  T)  x  U, 

(v)  E  x  n  =  0  on  (0,  T)  x  dQ, 

(vi)  E(0,x)  =  0,  H(0,x)  =  0  in  U. 


(9) 


(10) 


We  will  henceforth  drop  the 'symbol  on  D.  We  note  that  equation  (10,  iii)  follows  from 
the  continuity  equation  ^  y  .  j  =  o,  the  assumption  that  p(0)  =  0,  and  the  assumption 
that  V  •  Js  =  0  (in  the  sense  of  distributions).  The  modified  constitutive  laws  (2)  after 
substitution  of  Ohm’s  law  and  (9)  are 


(i)  D(f,x)  =  e0eoo(x)E(t,x)  +  Pr  +  J  <r(x)E(s,x)  ds, 
,  (ii)  B  =  ^0H, 

with  Pr  the  solution  of 

tPr  +  Pr  =  e0(es  -  e^E. 


(11) 


(12) 
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3  An  anisotropic  perfectly  matched  layer  absorbing  medium 


In  order  to  obtain  a  finite  computational  domain  we  surround  the  region  of  interest  with 
perfectly  matched  absorbing  layers  (PML’s).  These  layers  have  been  demonstrated  to  have 
many  orders  of  magnitude  better  absorbing  capabilities  than  many  of  the  approximate 
absorbing  boundary  conditions  [Taf98] .  In  this  section  we  will  first  construct  PML’s  that  can 
be  matched  to  free  space  and  then  we  will  discuss  how  this  construction  can  be  generalized 
to  the  case  of  a  Debye  medium. 

As  the  source  current  Js  is  zero  inside  the  PML  regions  we  will  neglect  Js  in  this  section. 
We  model  the  PML’s  as  lossy  regions  by  adding  artificial  loss  terms  to  Maxwell’s  equations. 
This  approach  is  based  on  using  anisotropic  material  properties  to  describe  the  absorbing 
layer  [SKLL95,  Ged96,  Rap95].  Starting  with  Maxwell’s  equations  in  the  time  domain  we 
will  derive  a  PML  model  in  the  frequency  domain  and  then  transform  the  corresponding 
equations  to  the  time  domain  by  taking  the  inverse  Fourier  transforms  of  the  frequency 
domain  equations.  We  define  the  modified  permeability  and  permittivity  tensors 


-  M  + 


Wm] 


e  =  e 


M. 

jw 


Using  the  definitions  (13)  we  define  two  new  constitutive  laws 


(13) 


Bnew  =  [/u]H, 

„  Dnew  =  [e]E, 


(14) 


where  for  every  field  vector  V,  V  denotes  its  Fourier  transform.  From  the  constitutive  laws 
(14)  we  can  write  Maxwell’s  equations  in  time-harmonic  form  as 


f  jcjDnew  —  V xH, 
j  CJ  B  new  —  V  XE, 


(15) 


v  •  Dnew  =  o, 

.  V  •  Bnew  =  0. 


The  split  field  PML  introduced  by  Berenger  [Ber94]  is  a  hypothetical  medium  based 
on  a  mathematical  model.  In  [MP95]  Mittra  and  Pekel  showed  that  Berenger’s  PML  was 
equivalent  to  Maxwell’s  equations  with  a  diagonally  anisotropic  tensor  appearing  in  the 
constitutive  relations  for  D  and  B.  For  a  single  interface,  the  anisotropic  medium  is  uniax¬ 
ial  and  is  composed  of  both  the  electric  and  magnetic  permittivity  tensors.  This  uniaxial 
PML  (UPML)  formulation  performs  as  well  as  the  original  split  field  PML  while  avoiding 
the  nonphysical  field  splitting.  As  will  be  shown  below,  by  properly  defining  a  general 
constitutive  tensor  [S'] ,  we  can  use  the  UPML  in  the  interior  working  volume  as  well  as  the 
absorbing  layer.  This  tensor  provides  a  lossless  isotropic  medium  in  the  primary  computa¬ 
tion  zone,  and  individual  UPML  absorbers  adjacent  to  the  outer  lattice  boundary  planes  for 
attenuation  of  spurious  wave  reflections.  The  fields  excited  within  the  UPML  are  also  plane 


wave  in  nature  and  satisfy  Maxwell’s  curl  equations.  The  derivation  of  the  PML  properties 
of  the  tensor  constitutive  laws  is  also  done  directly  by  Sacks,  et  al.,  in  [SKLL95]  and  by 
Gedney  in  [Ged96].  We  follow  the  derivation  by  Sacks,  et  al.,  here.  We  begin  by  considering 
planar  electromagnetic  waves  in  free  space  incident  upon  a  PML  half  space.  Starting  with 
the  impedance  matching  assumption,  i.e. ,  the  impedance  of  the  layer  must  match  that  of 
free  space:  £q  /io  =  [e] — 1  [/i]  we  have 

—  =  —  =  [S]  =  diag{oi,  ai,  a3}.  (16) 

eo  M  o 

Hence,  the  constitutive  parameters  inside  the  PML  layer  are  [e]  =  eo[*S']  and  [p]  =  Mo  [S']) 
where  [S']  is  a  diagonal  tensor. 

Using  (14)  and  (16),  we  thus  find  the  constitutive  laws  for  the  perfectly  matched  layer 
are 

Dnew  =  e0[S]E, 

(17) 

„  Bnew  —  Mo[S]H, 

where  the  tensor  [S]  is 

I"  a~1  0  O' 

[S]  =  0  a  0  .  (18) 

0  0  a 

The  perfectly  matched  layer  is  therefore  characterized  by  the  single  complex  number  a. 

Remark  1  To  match  the  anisotropic  PML  along  a  planar  boundary  to  a  dispersive  isotropic 
half-space,  we  will  need  to  modify  the  above  construction  for  the  phasor-domain  constitutive 
relation 

D  =  eoe^HE,  (19) 

where  e*(a>)  is  the  frequency  dependent  relative  permittivity  of  the  Debye  medium,  which 
will  be  defined  in  a  standard  manner  in  (31)  below. 

3.1  Implementation  of  the  uniaxial  PML  in  three  dimensions 

To  apply  the  perfectly  matched  layer  to  electromagnetic  computations,  we  replace  the  half 
infinite  layer  by  a  layer  of  finite  depth  backed  with  a  more  conventional  boundary  condition, 
such  as  a  perfect  electric  conductor  (PEC).  This  truncation  of  the  layer  will  lead  to  reflec¬ 
tions  generated  at  the  PEC  surface,  which  can  propagate  back  through  the  layer  to  re-enter 
the  computational  region.  In  this  case,  the  reflection  coefficient  R,  is  now  a  function  of  the 
angle  of  incidence  6,  the  depth  of  the  PML  S,  as  well  as  the  parameter  a  in  (18).  Thus,  this 
parameter  a  for  the  PML  is  chosen  in  order  for  the  attenuation  of  waves  in  the  PML  to  be 
sufficient  so  that  the  waves  striking  the  PEC  surface  are  negligible  in  magnitude.  Perfectly 
matched  layers  are  then  placed  near  each  edge  (or  face  in  the  3D  case)  of  the  computational 
domain  where  a  non-reflecting  condition  is  desired.  This  leads  to  overlapping  PML  regions 
in  the  corners  of  the  domain.  As  shown  in  [SKLL95] ,  the  correct  form  of  the  tensor  which 
appears  in  the  constitutive  laws  for  these  regions  is  the  product 

[5]  =  [sysysy  (20) 
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PML  PEC  WALL 


Figure  2:  PML  layers  surrounding  the  domain  of  interest.  In  the  corner  regions  of  the  PML, 
both  ax  and  az  are  positive  and  the  tensor  [S']  is  the  product  [Sj^fS]^.  In  the  remaining 
regions  only  one  of  az  (left  and  right  PML’s)  or  ax  (top  and  bottom  PML’s)  are  nonzero 
and  positive.  The  tensor  [S],  is  thus  either  [Sja,  or  [S]2  respectively.  The  PML  is  truncated 
by  a  perfect  electric  conductor  (PEC). 


where  component  [S]a  in  the  product  (20)  governs  attenuation  in  the  a  direction,  for  a  = 
x,  y ,  z  (see  Figure  2).  All  three  of  the  component  tensors  in  (20)  are  diagonal  and  have  the 
forms 

s"1  0  0  sy  0  0  sz  0  0 

[S]a=  0  s*  0  ;[S]y=  0  s~ 1  0  ;  [S]2  =  0  s2  0  .  (21) 

0  0  sx  0  0  sy  0  0  sj1 

In  the  above  sx,sy,  sz  are  analogous  to  the  complex  valued  parameter  a  encountered  in  the 
Section  3  analysis  of  the  single  PML  layer. 

When  designing  PML’s  for  implementation,  it  is  important  to  choose  the  parameters 
sa  so  that  the  resulting  frequency  domain  equations  can  be  easily  converted  back  into  the 
time  domain.  The  simplest  of  these  [Ged96]  which  we  employ  here,  is 

sa  =  1  H - — ,  where  aa  >  0,  a  =  x,y,z.  (22) 

J^eo 

The  PML  interface  represents  a  discontinuity  in  the  conductivities  aa.  To  reduce  the 
numerical  reflections  caused  by  these  discontinuous  conductivities,  the  aa  are  chosen  to  be 
functions  of  the  variable  a  (e.g.,  ax  is  taken  to  be  a  function  of  x  in  the  [5]^  component 
of  the  PML  tensor).  A  choice  of  these  functions  so  that  aa  =  0,  i.e.,  sa  =  1,  at  the 
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interface  yields  the  PML  a  continuous  extension  of  the  medium  being  matched  and  reduces 
numerical  reflections  at  the  interface.  If  one  increases  the  value  of  aa  with  depth  in  the 
layer,  one  obtains  greater  overall  attenuation  while  minimizing  the  numerical  reflections. 
Gedney  [Ged96]  suggests  a  conductivity  profile 


aa(a) 


Copt  |  oc  Oi  o 

5™ 


a  =  x,y,  z. 


(23) 


where  5  is  the  depth  of  the  layer,  a  =  ao  is  the  interface  between  the  PML  and  the 
computational  domain,  and  m  is  the  order  of  the  polynomial  variation.  Gedney  remarks 
that  values  of  m  between  3  and  4  are  believed  to  be  optimal.  Empirical  testing  suggests 
that,  for  a  broad  range  of  problems,  an  optimal  value  of  crmax  is  given  by 


m  +  1 

CTopt  ~  150ttAq^/M 


(24) 


where  Aa  is  the  space  increment  in  the  a  direction  and  er  is  the  relative  permittivity  of 
the  material  being  modelled.  In  the  case  of  free  space  er  =  1.  Since  we  desire  to  match  the 
PML  to  both  free  space  as  well  as  the  Debye  medium,  we  will  use  er  to  be  the  average  value 

er  =  ^(l  +  eoo)-  (25) 


3.2  Reduction  to  two  dimensions 


From  the  time-harmonic  Maxwell’s  curl  equations  in  the  PML  (15)  and  (17),  Faraday’s  and 
Amp  ere- Maxwell’s  laws  can  be  written  in  the  most  general  form  as 


ju;/iio[S']H  =  —  VxE,  (Faraday’s  Law) 

MS]  D  =  VxH  —  Js,  (Ampere-Maxwell’s  Law) 


(26) 
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In  (26),  [5]  is  the  diagonal  tensor  defined  via  (20),  (21)  and  (22).  In  the  presence  of  this 
diagonal  tensor,  a  plane  wave  is  purely  transmitted  into  the  uniaxial  medium.  The  tensor  [5] 
is  no  longer  uniaxial  by  strict  definition,  but  rather  is  anisotropic.  However,  the  anisotropic 
PML  is  still  referenced  as  uniaxial  since  it  is  uniaxial  in  the  non-overlapping  PML  regions. 
Let  H  =  X  x  Z  denote  the  computational  domain  along  with  the  absorbing  layers.  We 
partition  the  interval  X  into  disjoint  closed  intervals  as  X  =  A'Bpmi  U  Xq  U  ATpmi  and  the 
interval  Z  into  disjoint  closed  intervals  as  Z\jVm\  UZaU  Z^>  U  ZRpmi ,  as  seen  in  Figure  3.  The 
computational  domain  of  interest  is  the  region  Xq  x  {ZaU  Z\ 3},  where  Xq  x  Za  denotes  air 
and  Xq  x  Zjq  denotes  the  Debye  medium.  The  buffer  region  H  \  (Xq  x  {Za  U  Zjq})  contains 
the  absorbing  layers  (PML’s).  The  PML’s  are  backed  by  a  perfect  conductor  where  the 
boundary  condition  n  x  E  =  0  is  used.  To  obtain  a  two-dimensional  model,  we  make  the 
assumption  that  the  signal  and  the  dielectric  parameters  are  independent  of  the  y  variable. 
Also,  an  alternating  current  along  the  x  direction,  as  in  (45)  below,  then  produces  an  electric 
field  that  is  uniform  in  y  with  nontrivial  components  Ex  and  Ez  depending  on  (t,  x,  z )  which, 
when  propagated  in  the  xz  plane  result  in  oblique  incident  waves  on  the  dielectric  surface 
in  the  xy  plane  at  z  =  z\.  With  these  assumptions,  Maxwell’s  equations  reduces  to  a  two 
dimensional  system  for  solutions  involving  the  Ex  and  Ez  components  for  the  electric  field 
and  the  component  Hy  for  the  magnetic  field,  which  is  referred  to  as  the  TEy  mode.  In 
this  case,  we  have  oy  =  0  and  sy  =  1  in  the  UPML.  From  (11,  i)  and  (7)  expressed  in  the 
frequency  domain  we  have  the  constitutive  relation 


D  =  e0 


€s  ^OO 

1  +j  WT 


+  — ) 
jwe0y 


E. 


(27) 


Re-scaling  the  electric  field  as 

E=.f^E,  (28) 

V  IE 

we  can  write  the  tinre-harnronic  frequency  domain  Maxwell’s  equations  (26)  in  the  uniaxial 
medium  as 


(i) 

(ii) 

(hi) 


jw 

jw 


1  + 


1  + 


1  + 


<?z(z) 

jwe0 

<?x{x) 

jwe0 

ctx(x) 

jwe0 


1  + 


1  + 


1  + 


°x(x) 

jwe0 

<?z(z) 

j^e0 

<?z(z) 

j^e0 


-co( 


dHy 

dz 


+  Js  •  i), 


-1 


Dz  =  c0- 


dHy 

dx 


Hi, 


(29) 


In  the  above  i  is  the  unit  vector  in  the  x  direction,  co  =  3.0  x  108  m/s  is  the  speed  of  light 
in  vacuum  and  D  =  ( DX,DZ)T  is  defined  as 


with 


D  =  e*  E, 


^OO 


^00  ^  G 

1+jcur  jtue0' 


(30) 

(31) 
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To  avoid  a  computationally  intensive  implementation,  we  define  suitable  constitutive  rela¬ 
tionships  that  facilitate  the  decoupling  of  the  frequency  dependent  terms  [Taf98].  To  this 
end,  we  introduce  the  fields 

h*  —  s-1  f) 

^ X  —  °X  ^Xt 


D*  —  o-i n 
—  bz  ^zi 


H*y  =  SZHy. 


(32) 


Again,  dropping  the'symbol  on  D  and  transforming  the  frequency  domain  equations  to  the 
time  domain,  we  obtain  the  system 

F)  H*  i  u*  —  r  ( dEx 

B,H‘  +  ~To  ‘  ~  ^  (' 8^  “  8i" 

a,Hy  +  ^lHt  =  gtH^ 


eo 

(iii)  dtD*  4 - —D*  =  -co(—^-  +  Js  •  i) 

eo  oz 

(iv)  b,d,  -  a,Lr  +  , 

(v)  dtDt  +  E*lD:  =  CaaJk< 

eo  ox 

(vi)  dtDz  =  8tD*z  +  ^-D*z, 

eo 


(33) 


with 


D(f)  =  eooE(t)  +  P(t)  +  C(t),  (34) 

where  D  =  (Dx,  DZ)T ,  and  inside  of  the  dielectric  the  polarization  P  =  (PX,PZ)T  satisfies 

tP  +  P  =  (es  —  eoo)E,  (35) 

while  the  conductivity  term  C  =  (Cx,Cy)T  satisfies 

C  =  (a/e0)E  .  (36) 

inside  of  the  dielectric.  Outside  the  dielectric  we  have  P  =  C  =  0.  We  will  assume  zero 
initial  conditions  for  all  the  fields. 


3.3  Pressure  dependence  of  polarization 

We  consider  a  pressure  dependent  Debye  model  for  orientational  polarization  first  developed 
in  [ABR02] .  We  assume  that  the  material  dependent  parameters  in  the  differential  equation 
(35)  for  the  Debye  medium  depend  on  pressure  p,  i.e., 

r(p)P  +  P  =  (es{p)  -  eoo(p))E.  (37) 
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We  suppose  as  a  first  approximation  that  each  of  the  pressure  dependent  parameters  can 
be  represented  as  a  mean  value  plus  a  perturbation  that  is  proportional  to  the  pressure 

(i)  r(p)  =  r0  +  Krp, 

(ii)  es(p)  =  eSj0  +  Ksp,  (38) 

(hi)  Coo  (p)  =  £oo,0  +  KooP- 

This  can  be  considered  as  a  truncation  after  the  first-order  terms  of  power  series  expansions 
of  these  functions  of  p.  Thus  (35)  is  written  as 

(to  +  KTp) P  +  P  =  ((es,o  -  eoo,o)  +  («s  -  Koc)p)E.  (39) 

This  model  features  modulation  of  the  material  polarization  by  the  pressure  wave  and  thus 
the  behavior  of  the  electromagnetic  pulse.  We  will  assume  that  the  acoustic  pressure  wave 
is  specified  a  priori.  The  pressure  p  will  be  assumed  to  have  the  form 

(Z  —  Z2  \ 

up[t+— — ]J,  (40) 

where  Z2  G  Zd  with  z\  <  zi-  The  terms  up,  \p  and  cp  denote  the  acoustic  frequency, 
wavelength  and  speed  respectively.  We  note  that  outside  the  pressure  region  ( Z2,Z2  +  Ap), 
the  pressure  is  zero  and 

T  —  T0 ,  es  —  6s?o )  ^oo  =  ^oo,0  • 

From  (34) 

(i)  Dx(t)  =  (600^  +  KooP)Ex{t)  +  Px(t)  +  Cx(t) , 

(ii)  Dz(t)  =  (eoojO  +  Koop)Ez(i)  +  Pz(t)  +  Cz(t), 

with 

(1)  (ho  T  HtP)Px  T  Px  —  ((^s,0  ^(xjjO)  T  (ks  Koo'jp) Ex, 

(ii)  (t0  +  KTp)Pz  +  Pz  =  ((es,o  -  eoo>0)  +  (ks  -  n00)p)Ez, 

and 

f  (i)  Cx  =  (<r/e0 )EX, 

{  (ii)  Cz  =  (a/e0)Ez. 


(41) 

(42) 

(43) 

(44) 


3.4  Source  term 

The  source  term  Js  will  model  an  infinite  (in  the  y  direction)  antenna  strip,  finite  in  the  x 
direction  (between  x\  and  X2),  lying  in  the  z  =  zc  plane  in  free  space,  and  we  assume  the 
signal  is  polarized  with  oscillations  in  the  xz  plane  only,  with  uniformity  in  the  y  direction 
(see  again  Figure  1).  We  therefore  take 


Js(f,  x,  z)  =  I(Xl:X2)(x)5(z  -  zc)  sin(cuc(t  -  3f0))  exp 


t  -  3 10'  2\ 

to  \  J 


to 


2 

-  SGC,  UJC 

UJd  Wc 


27t fc  rad/sec. 


(45) 
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Here  I/XI  X2\  is  the  indicator  function  on  the  interval  {x\,X2)  £  Xq,  5(z  —  zc)  is  the  Dirac 
measure  centered  at  z  =  zc,  uc  is  the  central  angular  frequency  of  the  input  signal  and  cud 
is  the  decay  frequency.  The  Fourier  spectrum  of  this  pulse  has  even  symmetry  about  the 
central  frequency  fc.  The  pulse  is  centered  at  time  step  3to  and  has  a  1/e  characteristic  decay 
of  to  time-steps  [Taf95] .  We  will  specify  the  values  of  all  the  different  quantities  involved  in 
(45)  when  we  discuss  particular  problems  below.  The  source  (45)  radiates  numerical  waves 
having  time  waveforms  corresponding  to  the  source  function.  The  radiated  numerical  wave 
will  propagate  to  the  Debye  medium  and  undergo  partial  transmission  and  partial  reflection, 
and  the  solution  can  be  continued  until  all  transients  effectively  decay. 


3.5  Space  and  time  discretization 


We  will  use  the  FDTD  algorithm  [Taf95]  to  discretize  Maxwell’s  equations.  This  algorithm 
uses  forward  differences  to  approximate  the  time  and  spatial  derivatives.  The  FDTD  tech¬ 
nique  rigorously  enforces  the  vector  field  boundary  conditions  at  interfaces  of  dissimilar 
media  at  the  time  scale  of  a  small  fraction  of  the  impinging  pulse  width  or  carrier  period. 
This  approach  is  very  general  and  permits  accurate  modelling  of  a  broad  variety  of  materials 
ranging  from  living  human  tissue  to  radar  absorbers  to  optical  glass. 

The  system  of  equations  (33)  can  be  discretized  on  the  standard  Yee  lattice.  The  domain 
D  is  divided  into  square  cells  of  length  h  =  Ax  =  A z,  h  being  the  spatial  increment.  The 
electric  field  degrees  of  freedom  are  at  the  midpoints  of  edges  of  the  squares  and  the  magnetic 
degrees  of  freedom  are  at  the  centers  of  cells.  The  time  interval  over  which  time  stepping 
is  done  is  divided  into  subintervals  of  equal  size  using  the  time  increment  At.  We  perform 
normal  leapfrogging  in  time  and  the  loss  terms  are  averaged  in  time.  This  leads  to  an 
explicit  time  stepping  scheme  (see  [BB04b]  for  details)  .  In  this  case,  a  stability  condition 
(CFL)  has  to  be  satisfied  [Taf95]  in  order  to  obtain  a  well  posed  computational  scheme.  The 
time  step  At  and  the  spatial  increment  h  for  the  FDTD  scheme  in  non-dispersive  dielectics 
are  related  by  the  condition 


n cn 


coAf  1 

h  <  ^ 2 


(46) 


The  number  rj cn  is  called  the  Courant  number.  We  fix  the  value  of  r] cn  =  1/2.  In  [Pet94], 
the  author  established  that  several  extensions  of  the  FDTD  for  Debye  media  satisfy  the 
stability  restriction  of  the  standard  FDTD  scheme  in  nondispersive  dielectrics.  He  has  also 
shown  that  the  extensions  do  not  preserve  the  non-dissipative  character  of  the  standard 
FDTD  scheme.  The  extended  difference  systems  are  more  dispersive  than  the  standard 
FDTD,  and  their  accuracy  depends  strongly  on  how  well  the  chosen  timestep  resolves  the 
shortest  timescale  in  the  problem  regardless  of  whether  it  is  the  incident  pulse  timescale, 
the  medium  relaxation  timescale,  or  the  medium  resonance  timescale. 


4  The  forward  problem  for  a  Debye  medium:  First  test  case 

As  discussed  in  Section  3.3  the  pressure  dependent  parameters  es,  e^,  and  r  are  represented 
as  a  mean  value  plus  a  perturbation  that  is  proportional  to  the  pressure  in  equations  (38) 
(as  discussed  in  [ABR02]).  In  our  first  test  case,  we  present  numerical  simulations  for  a 
Debye  medium  that  is  similar  to  water.  For  signal  bandwidths  in  the  microwave  regime 
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the  dispersive  properties  of  pure  water  are  usually  modelled  by  a  Debye  equation  having  a 
single  molecular  relaxation  term.  The  mean  values  e*  0,  0,  and  Tq  for  this  Debye  model 

[APM89]  are  given  by 


e*  0  =80.1  (relative  static  permittivity), 

eoo  o  =  5.5  (relative  high  frequency  permittivity), 

Tq  =  8.1  x  10'12  seconds, 

a*  =  1  X  10”5  mhos/meter, 


(47) 


where  the  *  superscript  will  denote  the  true  values  of  all  the  corresponding  dielectric  pa¬ 
rameters.  Since  we  do  not  yet  have  experimental  data  to  determine  the  actual  values,  we 
instead  choose  trial  values  for  the  coefficients  of  pressure,  n*,  k ^  and  n*  to  use  in  generating 
simulated  data.  As  a  first  approximation  we  take  these  trial  values  to  be  a  fraction  of  the 


mean  values  of  the  dielectric  parameters,  e*  0,  0,  and  Tq  ,  respectively.  Thus, 

we  choose 

< 

=  0.6e*0  =  48.06, 

K*oc 

=  0.8e^0  =  4.4, 

(48) 

< 

=  0.05tq  =  4.05  x  10~13. 

We  hope  to  determine  the  values  of  these  pressure  coefficients  from  appropriately  designed 
experiments  in  the  near  future.  We  have  assumed  that  the  effect  of  the  electric  field  on  the 
acoustic  pressure  is  negligible.  We  are  interested  in  how  and  to  what  extent  the  acoustic 
pressure  can  change  the  reflected  electric  wave  from  the  interface  at  z  =  zi.  Moreover, 
the  relaxation  parameter  is  a  characteristic  of  the  material  in  [21,22],  which  affects  the 
transmission  of  electric  waves  through  this  region.  Consequently,  we  will  examine  in  more 
detail  how  tq  can  affect  the  electric  field.  The  electromagnetic  input  source  has  the  form 


J s(t,  x,  z)  =  I(XljX2)S(z)  sin (uc(t  -  3 10))  exp  \  - 

1 


t  —  3fo 
to 


(49) 


to  — 


2ir  x  109 


sec,  ujc  =  67T  x  lCr  rad/sec,  fc  =  3.0  x  10J  Hz. 


The  Fourier  spectrum  of  this  pulse  has  even  symmetry  about  3.0  GHz.  The  computational 
domain  is  defined  as  follows.  We  take  Xq  =  (0,0.1),  Za  =  (0,0.15)  and  Zd  =  (0.15,0.2). 
The  number  of  nodes  along  the  2-axis  is  taken  to  be  320  and  the  number  of  nodes  along  the 
x-axis  is  taken  to  be  160.  The  spatial  step  size  in  both  the  x  and  2  directions  is  Ax  =  A 2  = 
h  =  0.1/160.  From  the  CFL  condition  (46)  with  the  Courant  number  77 cn  =  1/2  we  obtain 
the  time  increment  to  be  At  ~  1.0417  pico  seconds.  The  central  frequency  of  the  input 
source  as  described  in  (49)  is  3.0  GHz  and  based  on  the  speed  of  light  in  air,  co  =  3x  108  m/s, 
we  calculate  the  corresponding  central  wavelength  to  be  Ac  =  {2ttcq)/u  =  0.1  meters.  The 
antenna  is  half  a  wavelength  long  and  is  placed  at  (xi,X2)  x  zc,  with  zc  =  0,  x\  =  0.025 
and  X2  =  0.075.  We  use  PML  layers  that  are  half  a  wavelength  thick  on  all  fours  sides 
of  the  computational  domain  as  shown  in  Figure  3.  The  reflections  of  the  electromagnetic 
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pulse  at  the  air-Debye  interface  and  from  the  acoustic  pressure  wave  are  recorded  at  the 
center  of  the  antenna  (xc,  zc),  with  xc  =  0.05,  at  every  time  step.  This  data  will  be  used  as 
observations  for  the  parameter  identification  problem  to  be  presented  later.  The  component 
of  the  electric  field  that  is  of  interest  here  is  the  Ex  component.  Thus  our  data  is  the  set 
E(q*)  =  {Ex(nAt,  xc,  zc;  q*)}^,  where  q*  =  (e*  0,  e^  0,  r0*,  a*,  k*s,  k*)t.  The  windowed 


acoustic  pressure  wave  as  defined  in  (40)  has  the  parameter  values  uip  =  6.07T  x  105  rad/sec, 
cp  =  1500  nr/s,  and  thus,  Xp  =  0.005  nr.  The  location  of  the  pressure  region  is  in  the  interval 
(z2,Z2  +  \p)  =  (0.175,0.18). 

In  Figure  4  (left)  we  plot  the  electromagnetic  source  that  is  used  in  the  simulations  for 
the  Debye  model  (47).  We  plot  the  power  spectral  density  of  the  source  (49)  in  Figure  4 
(right).  The  power  spectral  density  |y|2  of  the  vector  Js  is  defined  to  be 


Z  =  FFT(JS), 
T|2  =  Z-Z, 


(50) 


where  FFT(JS)  is  the  fast  fourier  transform  of  the  vector  Js.  As  seen  in  Figure  4  (right) 
the  power  spectral  density  is  symmetric  about  3.0  GHz. 

In  Figure  5  we  plot  the  Ex  field  magnitude  at  the  center  of  the  antenna  versus  time, 
which  shows  the  electromagnetic  source,  the  reflection  off  the  air-Debye  interface  and  the 
reflection  from  the  region  containing  the  acoustic  pressure  wave.  As  seen  in  these  plots  the 
amplitude  of  the  reflection  from  the  acoustic  pressure  is  many  orders  of  magnitude  smaller 
than  that  of  the  initial  electromagnetic  source  as  well  as  that  of  the  reflection  from  the 
air-Debye  interface.  Thus,  it  is  an  interesting  question  as  to  whether  the  reflections  from 
the  acoustic  pressure  region  can  be  used  for  identification  of  parameters. 

In  Figure  6  we  plot  the  Ex  field  magnitude  at  different  times  in  the  plane  containing 
the  center  of  the  antenna  versus  depth  along  the  z  axis.  The  top  left  figure  shows  the 
electromagnetic  wave  penetrating  the  Debye  medium.  In  the  top  right  figure  we  see  the 
reflection  of  the  electromagnetic  wave  from  the  Debye  medium  moving  towards  the  antenna 
and  the  Brillouin  precursor  propagating  in  the  Debye  medium.  In  the  bottom  left  figure  we 
observe  the  reflection  from  the  region  containing  the  acoustic  pressure  and  the  transmitted 
part  of  the  electromagnetic  source  travelling  into  the  absorbing  layer.  The  reflection  from 
the  acoustic  pressure  crosses  the  Debye  medium  in  the  bottom  right  figure. 


4.1  Sensitivity  analysis 


In  this  section  we  examine  the  system  dynamics  as  the  parameters  vary.  In  particular,  we 
are  interested  in  the  changes  produced  by  the  coefficients  of  pressure  in  the  polarization, 
namely  the  parameters  ns,  Koo  and  kt. 

We  consider  equation  (31)  to  argue  that  the  parameter  ks  appears  to  be  the  most 
influential  in  the  wave  interaction,  whereas  the  pressure  coefficients  Kqq  and  nT  do  not 
seem  to  be  as  important.  In  our  simulations  the  outgoing  and  reflected  radiation  will  be 
dominated  by  frequencies  near  the  center  frequency  3.0  GHz.  Thus,  e*  will  be  dominated  by 
frequencies  near  3.0  GHz.  In  this  problem  cutq  ~  0(  10~2),  and  eo  =  8  x  10-12.  Rewriting 
(31)  as 

er  —  TT~- - ^  (  TT^ —  )  e°°  +  : - ’ 

1+jtUT  \l+]UJTj  jcueo 
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x  i  o3  The  power  spectral  density  of  the  Source 


Figure  4:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (47)  and  (48). 
(Left)  The  source  as  defined  in  (49).  (Right)  Fourier  transform  of  the  source.  The  transform 
is  centered  around  the  central  angular  frequency  of  luc  =  6tt  x  109  rad/sec. 


Figure  5:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (47)  and  (48). 
Time  plot  of  the  Ex  component  of  the  electric  field.  This  is  the  ’’data”  that  is  received  at  the 
center  of  the  antenna.  (Left)  The  original  signal  (source),  the  reflection  off  the  Air-Debye 
interface  and  the  reflection  due  to  the  pressure  wave.  (Right)  The  (scale  enlarged)  reflection 
off  the  Air-Debye  interface  and  the  reflection  due  to  the  pressure  wave. 
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t=  1.04  x  10  9  sec 


t  =  1 .46  x  1 0  9  sec 


t=  2.08  x  10  9  sec 


t  =  2.5x  1 0  9  sec 


Figure  6:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (47)  and 
(48).  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z  (meters). 
(Top,  left)  The  input  signal  gets  partially  reflected  off  the  Debye  medium  and  partially 
transmitted  into  the  medium.  (Top,  right)  The  reflection  off  the  air-Debye  interface  is  seen 
moving  towards  the  antenna,  while  the  transmitted  part  of  the  source  is  seen  propagating 
in  the  Debye  medium  towards  the  region  containing  the  pressure  wave.  (Bottom,  left)  The 
interaction  of  the  electromagnetic  wave  with  the  pressure  wave  causes  the  electromagnetic 
wave  to  partially  reflect  and  partially  transmit.  (Bottom, right)  The  reflection  from  the 
acoustic  pressure  wave  impinges  on  the  air-Debye  interface,  while  the  wave  transmitted  into 
the  pressure  region  moves  into  the  right  absorbing  layer. 
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we  consider  each  of  the  three  terms  in  (51)  separately  to  determine  their  magnitude.  The 
magnitudes  of  the  corresponding  true  values  of  these  three  terms  (neglecting  the  acoustic 
pressure  terms)  are  approximately  given  as 


fcs,0 


1  +j  WT0* 

1  +j^T0  / 

ex’ 


COO,0 


J^e0 


O(e*0) 

«  0(10), 

0(  10-2e^0) 

«  e>(io-2 

e>(ioV) 

«  e>(io-3 

(52) 


Thus  we  see  that  e*  will  be  most  sensitive  to  the  static  permittivity  es,o  and  the  effects  of 
600,0  and  ro  will  not  be  as  pronounced.  Moreover,  this  implies  that  e*  will  be  more  sensitive 
to  the  pressure  coefficient  ns  than  to  the  coefficients  Koo  and  kt. 

We  are  also  interested  in  determining  the  effect  that  the  acoustic  speed  and  frequency 
have  on  the  amplitude  of  the  reflection  from  the  acoustic  pressure  region.  By  changing 
the  speed  or  the  frequency,  the  wavelength  Xp  changes  and  thus  the  size  of  the  interval 
(z2,Z2  +  Xp)  in  which  the  pressure  wave  is  generated.  In  [BB04b]  we  have  shown  that 
the  amplitude  of  the  acoustic  reflection  appears  to  decrease  as  the  acoustic  frequency  is 
increased  and  increases  at  the  acoustic  speed  increases.  We  note  here  that  our  windowed 
pressure  wave  contains  only  one  wavelength  of  the  sinusoid.  Thus  in  general  it  will  be 
difficult  to  predict  precisely  how  the  reflections  will  behave  as  the  speed  or  the  frequency  is 
changed. 


5  Parameter  estimation  and  statistical  inference:  The  in¬ 
verse  problem 

In  our  forward  simulations  the  signal  that  we  record  at  the  center  of  the  antenna  is  a 
set  of  measurements  of  Ex,  the  x  component  of  the  electric  field,  at  the  point  (xc,zc)  in 
our  computational  domain  at  discrete,  uniformly  spaced  intervals  of  time.  This  signal  has 
segments  consisting  of  the  input  source,  the  reflection  at  the  air-Debye  interface  and  another 
reflection  from  the  region  that  contains  the  pressure  wave  as  seen  in  Section  4.  The  reflected 
signal  is  a  function  of  the  various  dielectric  properties  of  the  Debye  medium;  namely,  the 
static  relative  permittivity  es  o,  the  infinite  frequency  permittivity  600,0,  the  relaxation  time 
ro,  the  conductivity  a,  and  the  three  pressure  coefficients  ks,  Koo  and  kt.  We  collect  all 
these  quantities  together  to  define  the  parameter  vector 

Q  —  [fs,0,  £00, 0,  7"0,  ^OO,  ^t]  •  (53) 

Thus,  the  signal  recorded  in  the  forward  simulation  is  a  function  of  q.  In  general,  such  sig¬ 
nals  are  usually  obtained  in  an  experimental  setting  conducted  in  a  laboratory  with  physical 
equipment  such  as  electric  pulse  generators  to  generate  electric  pulses,  piezoelectric  trans¬ 
ducers  that  generate  the  acoustic  pressure  waves  and  antennas/receivers  that  record  the 
electric  field  intensities.  Moreover,  such  signals  generated  in  an  experimental  setting  are 
usually  noisy,  with  noise  arising  through  the  equipment  that  is  used;  namely  the  resistor, 
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antenna/receiver  and  the  transducer,  in  our  case.  Implicit  in  modelling  such  a  collection  of 
measurements  is  the  assumption  that  there  exist  true  values  of  all  the  parameters  that  char¬ 
acterize  the  medium  to  be  interrogated.  We  will  denote  the  corresponding  true  parameter 
vector  by  q* . 

In  order  to  simulate  an  experimental  situation  (i.e. ,  generate  typical  “data”),  we  use 
these  true  values  (which  of  course  would  not  be  known  in  an  experimental  setting!)  in 
a  forward  simulation ,  and  record  the  signal  observed  at  the  antenna  as  described  above. 
We  then  add  noise  to  the  generated  signal  to  produce  a  noisy  signal  which  will  form  the 
observations  or  data.  We  use  simulated  data  in  the  inverse  problem  for  estimation  of  q*  to 
validate  the  methods  first  on  “data”  from  known  parameters  in  a  setting  with  known  noise. 
If  we  cannot  estimate  the  dielectric  parameters  from  data  that  is  constructed  via  numerical 
simulations  of  our  discrete  model,  then  the  reconstruction  of  these  parameters  from  actual 
experimental  data  usually  will  not  be  feasible. 

We  can  state  the  inverse  or  parameter  estimation  problem  that  we  will  attempt  to  solve 
as  follows:  Using  observations  or  data  (electric  field  intensities  containing  noise  collected 
at  the  center  of  the  antenna),  determine  an  estimate  q,  of  the  true  parameter  q*,  that 
belongs  to  an  admissible  set  Q,  so  that  the  solution  of  the  forward  problem  with  q  =  q  best 
describes  the  data  that  is  collected.  When  we  have  solved  this  deterministic  problem  (as  will 
be  detailed  below),  it  will  also  be  necessary  to  specify  reliability  of  our  estimates,  i.e.,  can 
we  specify  measures  of  uncertainty  related  to  the  estimate  q  of  q*?  It  is  important  to  note 
that  the  uncertainty  in  consideration  is  inherent  in  the  method  of  producing  the  estimates 
as  well  as  in  the  process  of  data  collection.  Such  measures  of  uncertainty  will  be  specified  by 
means  of  confidence  intervals.  These  intervals  are  a  probability  statement  about  the  inverse 
problem  and  the  procedure  which  is  used  to  construct  estimates  of  the  parameters.  Thus  we 
are  led  in  a  completely  natural  way  to  stochastic  or  probabilistic  aspects  of  estimates  from 
a  deterministic  problem  solved  with  deterministic  algorithms  [BB].  The  statistical  error 
analysis  that  we  will  present  here  is  based  on  standard  statistical  formulations  as  given  in 
[DG95]. 

We  first  consider  the  two  different  ways  in  which  we  can  add  noise  to  our  signal,  i.e.,  the 
values  of  the  x  component,  Ex,  of  the  electric  field  observed  at  the  center  of  the  antenna, 
(xc,  zc),  for  different  times  tk  =  kAt,  k  =  1, . . . ,  M.  We  represent  this  signal  as  the  vector 

E(q*)  =  {E*k}l li  =  {Ex(tuxc,  zc]  q*), . . . , Ex(tM,xc,  zc]  q*))T  ,  (54) 

where  q*  =  (e*  0,  0,  Tg ,  a*,  k*,  k*)t  is  the  true  parameter  vector. 

1.  Relative  random  noise  (RN)  :  In  this  case  the  amplitude  of  the  noise  that  is  added  is 
proportional  to  the  size  of  the  signal,  {E^}^fL1.  Our  simulated  data  samples  are 

0(  =  Eld  +  vr,%),  k  =  1 , . . .  M,  (55) 

where  fj  =  {%}^£1  are  independent  normally  distributed  random  variables  with  mean 
zero  and  variance  one,  i.e.,  rjk  ~  jV(0,  l),fc  =  1 , . . . ,  M.  We  express  the  relative 
magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  E(q*)  by  taking  two  times 
the  value  of  u  as  the  size  of  the  random  variable.  For  example,  v  =  0.005  corresponds 
to  1%  relative  noise  [BBL00].  This  noise  model  does  not  produce  constant  variance 
across  the  samples. 
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2.  Constant  variance  noise:  Since  constant  variance  is  most  conveniently  assumed  in 
standard  error  analysis,  we  will  consider  estimates  obtained  from  an  inverse  problem 
applied  to  simulated  data  which  contains  constant  variance  random  noise.  The  data 
that  we  try  to  fit  in  this  case  is 

Osk  =  El  +  f3r,l  k  =  1  .  M,  (56) 

where  as  in  (55),  77*.  ~  jV(0,  l),k  =  1  The  constant  (3  is  taken  to  be  the 

product  va,  where  a  is  a  scaling  factor  which  is  chosen  so  that  the  signal  to  noise 
ratio  (SNR)  of  the  data  defined  in  (56)  corresponds  to  that  of  the  data  defined  in  (55) 
with  noise  level  v.  This  justifies  comparing  results  obtained  with  the  two  different 
ways  of  adding  comparable  noise  to  our  signal.  This  will  be  addressed  in  more  detail 
in  Section  6.2. 

The  vector  Os  =  {Ol}¥  l  will  be  our  data  or  observations,  and  the  noise  fjs  =  {r]t}uL1  will 
be  referred  to  as  the  measurement  errors  corresponding  to  the  observations.  The  statistical 
error  analysis  that  we  now  develop  is  only  applicable  to  the  case  of  constant  variance  noise. 
Thus,  for  the  rest  of  this  section,  we  will  assume  that  our  observations  are  of  the  form 
(56).  The  sample  observations  Os,  which  are  in  general  obtained  from  experiment,  are  a 
realization  of  the  corresponding  random  variable  O  =  (O i,  O 2, . . . ,  Om)T  which  can  be  seen 
to  be  a  transformation  of  the  random  variable  fj.  Thus 

Ok  =  E*k  +  (3Vk,  k  =  l,...M,  (57) 

is  a  stochastic  process  and  has  a  multivariate  normal  distribution  with  mean  vector  E(q*) 
defined  in  (54)  and  covariance  matrix  (32Im ■  The  statistical  model 

0~A7M{E(q*),/32TM},  (58) 

is  a  formal  representation  of  the  population  of  all  possible  realizations  of  0 1,  O2,  ■  ■  ■ ,  Om  that 
can  be  observed.  When  we  collect  data,  we  are  observing  a  single  realization  Of, ... ,  OsM, 
i.e.,  a  sample.  Our  objective  then  is  to  estimate  the  true  value  of  the  parameter  q*  by 
collecting  data,  i.e.,  observing  a  single  realization  of  O  as  well  as  accounting  for  the  fact 
that  a  different  realization  will  produce  a  different  estimate  of  q*.  We  would  like  to  learn 
about  the  true  value  q*  (which  determines  the  nature  of  the  population)  from  a  sample, 
as  well  as  indicate  the  certainty  (or  uncertainty)  that  we  can  associate  to  our  knowledge  of 
this  parameter.  This  process  of  making  statements  about  a  population  of  interest  on  the 
basis  of  a  sample  from  the  population  is  called  statistical  inference. 

The  unknown  true  parameter  vector  q*  will  be  estimated  by  means  of  a  suitable  function 
q(O)  of  the  observations  O.  The  function  q(O)  is  a  random  variable  and  is  referred  to  as 
an  estimator.  When  evaluated  at  a  particular  realization,  Of ,  Of, . . . ,  OsM,  this  estimator 
yields  a  numerical  value  that  gives  information  about  the  true  value  of  the  parameter  q*. 
For  a  particular  realization  Os  of  the  random  variable  of  observations  O,  we  call  q(Os)  an 
estimate.  In  this  paper  we  will  consider  the  method  of  least  squares  to  obtain  estimates 
for  the  parameters  and  to  calculate  confidence  intervals  for  the  estimates  by  linearizing 
around  the  estimate  q(Os).  In  this  case  q(O)  will  be  the  least  squares  estimator  qoLs(O)- 
Associated  with  the  random  variable  qoLs(O)  is  the  probability  space  Q  =  where 
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Q  is  the  sample  set  of  all  admissible  parameters  q,  B  is  the  cr-algebra  of  events  and  m  is  the 
probability  measure  (or  distribution).  Thus  we  look  for  the  least  squares  estimator  qoLs(O) 
on  Q  such  that 


qoLs(O)  =  arg  rninJ(q), 
q  e  Q 

1  M 

J(q)  =  2  X!  I Ex{tk,xc  zc;  q)  -  ok |2 , 

A  k=  1 


(59) 


where  E(q)  =  (Ex(ti,xc,  zc ;  q),  Ex(t2,  xc,  zc ;  q), . . . ,  Ex(tM,xc,  zc ;  q))T  will  be  generated  by 
our  forward  simulation.  We  will  refer  to  the  vector  E(q)  as  our  simulations.  For  a  particular 
realization  Os  of  the  random  variable  O  of  observations,  this  process  will  yield  the  least 
squares  estimate  qoLs(Os)  £  Qi  where 

qoLs(Os)  =  arg  min  J5(q), 
qe<3 

1  m  (60) 

jS(q)  =  2  l^(4,a;c,^c;q)  -  osk \2 . 

z  k=  1 


5.1  Calculation  of  confidence  intervals  via  linearization 


In  general,  the  mappings  between  the  parameters  q  and  the  simulations  E  are  nonlinear. 
Let  Ek  =  Ex(tk,  xc,  zc',  q).  The  functions 


Ek{qi, . . .  ,qi)  =  Ex(tk,xc,zc-,q),  k=l,...,M,  (61) 

denote  real-valued  differentiable  functions  of  the  unknown  parameters  q  =  (qi, . . . ,  qi)T , 
where  1  <  l  <  7  depending  on  how  many  parameters  we  attempt  to  identify.  Let  q  = 
q*  +  T’q,  where  the  corrector  term  X>q  =  (Agi, . . .  Aqi)T  is  unknown.  Linearizing  the 
functions  in  (61)  around  the  true  parameter  q*  by  using  the  Taylor  expansion,  we  obtain 


Let  us  define 

and 


Ek(qi,  •••<?/)  =  Ek{q{  +  Aqu...,qi  +  A  qi) 

i 


A  q.j ,  k  =  1, ,  M. 


*\\T 


Vc  =  (Oi  -  ^(q*),  Em{  q*)) 


^(q*)  = 


/ 

dEi  1 

dEi 

dqi  |  q* 

dqi 

8Em  I 

9Em 

V 

dqi  1  q* 

dqi 

(62) 


(63) 


(64) 


) 


Thus,  in  order  to  estimate  the  unknown  corrector  terms  by  the  method  of  least  squares,  the 
estimator  Pols  is  determined  by 


1 


Pols  =  min  -(Pc  -  ;f(q*)Pq)J  (Pc  -  *(q*)Pq), 


vq  2 


(65) 
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which  is  the  linearized  version  of  (59).  This  involves  solving 


8  (peTPe  -  2PeTT’(q*)P0Ls  +VTLSX(q*)TX(cf)VoLs) 

<9£>ols 

and  the  corresponding  optimality  condition  yields  the  estimator 


=  0, 


(66) 


Adls  =  (df(q*)Tdf(q*))-1df(q*)TPe.  (67) 

We  note  that  the  expected  value  of  Pols  is  Exp(PoLs)  =  0)  anci  the  covariance  matrix  of 
Pols  and  hence  also  of  qoLS  is 

Cov(Pols)  =  (df(q*)Tdf(q*))-1df(q*)TCov(Pe)T’(q*)(T’(q*)TT’(q*))~1 

O  rr  ,  (68 

=  (3\X{c{)TX{cC)rl  =  Cov(qoLs), 


where  qoLS  =  q*  +  Pols>  and  f3  is  defined  in  (56).  Thus  the  least  squares  estimator  is  a 
random  variable  with  a  multivariate  normal  distribution  which  after  linearization  can  be 
approximately  represented  as 

qoLs(O)  ~  Mi  (q*,/32[T’T(q*)T’(q*)]-1)  .  (69) 

However,  when  data  is  generated  in  an  experimental  setting,  we  do  not  have  any  knowledge 
of  the  true  parameter  vector  q*  and  hence  we  cannot  calculate  (df(q*)TT’(q*))  .  Thus, 

we  further  approximate  our  least  squares  estimator  as  having  the  multivariate  normal  dis¬ 
tribution 

qoLs(O)  ~A fi  (qoLs(Os),^Ls[^T(qOLs(Os))df(qOLs(Os))]-1)  ,  (70) 

which  is  obtained  by  repeating  the  linearization  process  around  the  estimate  qoLs(Os) 
obtained  from  a  realization  Os  of  the  random  variable  O  of  observations.  The  value  of 
$ ols  chosen  as 


Pols  =  £  I Ex(*k,  «c;  qoLs(Os))  -  Osk |2 ,  (71) 

L  k=  1 

which  is  the  minimum  value  of  the  least  squares  objective  function  scaled  by  2/{M  —  l ). 

Whenever  an  estimate  based  on  data  is  reported,  it  should  be  accompanied  by  an  as¬ 
sessment  of  uncertainty  based  on  the  sampling  distribution.  To  this  end  we  will  construct 
confidence  intervals  for  all  the  estimates  that  we  will  provide.  The  approach  we  use  is  to 
look  at  the  standard  error  (SE)  for  each  of  the  l  components  of  the  estimate  qoLs(Os)> 
which  is  given  for  its  jth  component  by  the  jth  diagonal  term  of  the  covariance  matrix,  i.e., 


SE(c/0Ls,j)  =  \/var(g0LSj(Os)) 

V, -  (72) 

=  \/P2OLs((X^OLSmrX^OLSm))-%  j  = 
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We  construct  the  intervals 


CIj  =  (<7oLS,j(Os)  -  1-96  SE(qoLSj),  QOLS,j{Os)  +  1.96  SE(q0LS,j)) ,  j  =  1,  ■  •  •  J,  (73) 
for  which 

m  ({^OLS,j(Os)  -  1-96  SE(q0LS,j)  <  90LS,j(0)  <  <ZOLSj(Os)  +  1.96  SE(qOLs,j)})  =  0.95, 

j  =  1,-  •  •  ,h 

(74) 

where  m  is  the  probability  measure  for  the  probability  space  Q.  and  M  is  sufficiently  large 
that  one  can  use  a  Gaussian  A7(0, 1)  distribution  in  computing  confidence  intervals.  These 
intervals  CIj  are  called  confidence  intervals  with  confidence  level  0.95  or  the  95%  confidence 
interval. 

Remark  2  The  confidence  intervals  are  a  probability  statement  about  the  procedure  by 
which  an  estimate  is  constructed  from  a  sample  of  the  population.  They  are  not  a  prob¬ 
ability  statement  about  the  true  parameter  q*  which  is  fixed.  If  we  could  construct  the 
confidence  intervals,  from  our  estimation  procedure,  for  all  possible  data  samples  of  size  M, 
then  95%  of  such  intervals  would  cover  the  true  parameter  values  q*.  However,  for  a  par¬ 
ticular  data  sample  Os,  the  confidence  intervals  constructed  as  above  may  or  may  not  cover 
the  true  parameter  q* .  What  we  can  state  is  that  we  are  95%  confident  that  the  confidence 
intervals  constructed  by  our  estimation  procedure  would  cover  q* . 


6  Parameter  estimation:  First  test  problem 

We  now  attempt  to  estimate  the  dielectric  parameters  for  the  Debye  medium  defined  by 
(47)  and  (48).  In  particular  the  values  of  all  the  parameters  are 


f* 

es,0 

=  80.1 

c* 

"oo,0 

=  5.5 

T* 

To 

=  8.1  x  10‘ 

a* 

=  1  x  10“5 

k*s  =  48.06 

<c  =4-4 

k*  =  4.05  x  10“13. 


(75) 


We  will  refer  to  the  values  (75)  as  the  true  values  of  the  parameters.  We  solved  the  corre¬ 
sponding  least  squares  optimization  problem  using  two  different  methods;  the  Nelder-Mead 
algorithm,  which  is  a  simplex  based,  gradient  free  method,  and  the  Levenberg-Marquardt 
trust  region  method  that  uses  forward  finite  differences  to  calculate  the  gradient.  From  our 
analysis  in  Section  4.1,  we  found  that  our  problem  is  sensitive  to  only  three  of  the  seven 
parameters.  We  will  then  try  to  estimate  two  or  three  of  the  parameters  to  which  our 
problem  is  sensitive.  The  particular  implementation  of  the  Levenberg-Marquardt  and  the 
Nelder-Mead  algorithm  that  are  used  for  the  calculations  presented  in  this  paper  are  based 
on  the  corresponding  algorithms  given  in  [Kel99] 
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Parameter 

True  Values 

Final  Estimates 

e«,o 

80.1 

80.1111 

^oo,0 

5.5 

5.7461 

TO 

8.1  x  10~12 

8.1283  x  10~12 

a 

1.0  x  10”5 

1.1972  x  10-5 

Ks 

48.06 

48.1227 

Hoo 

4.4 

4.9148 

K,t 

4.05  x  10~13 

4.1935  x  10~13 

Table  1:  Final  estimates  for  all  seven  parameters 
for  water  using  Nelder-Mead.  The  initial  param¬ 
eter  simplex  has  values  that  are  5-10%  larger 
than  the  true  values. 


6.1  Simulation  results:  The  Nelder-Mead  method 


We  first  present  parameter  estimation  results  using  the  Nelder-Mead  simplex  based  algo¬ 
rithm  which  is  used  for  optimization  in  noisy  problems.  The  Nelder-Mead  simplex  algorithm 
keeps  a  simplex  S  of  approximations  to  an  optimal  point.  In  this  algorithm  the  vertices 
of  the  simplex  {q?.}^ij,  where  l  is  the  size  of  q,  are  sorted  according  to  the  least  squares 
objective  function  values 

Js(qi)  <  Js(q2)  <  •••  <  Js(q;+i)-  (76) 


The  point  qi  is  called  the  best  vertex  and  q;+i  the  worst.  The  algorithm  attempts  to 
replace  the  worst  vertex  qz+i  with  a  new  point  of  the  form 


q  {ynm)  —  (IT  ^nm)q  ^nmqi+b 
where  q  is  the  centroid  of  the  convex  hull  of  {qi}(=1 


i= 1 


(77) 


(78) 


The  value  of  vnm  is  selected  from  a  sequence  of  values  that  are  computed  by  performing 
several  different  operations  on  the  simplex.  In  general,  the  Nelder-Mead  algorithm  is  not 
guaranteed  to  converge,  even  for  smooth  problems.  The  failure  mode  is  stagnation  at  a 
nonoptimal  point.  For  further  details  we  refer  the  reader  to  [Kel99]. 

In  the  first  test  we  will  try  to  estimate  all  of  seven  parameters  in  q* .  The  initial  simplex 
is  chosen  to  have  values  that  are  5-10%  larger  than  the  true  parameter  values  given  in  (75). 
We  do  not  add  any  noise  to  our  data  in  this  first  test.  The  final  estimates  from  the  Nelder- 
Mead  algorithm  are  presented  in  Table  1.  The  Nelder-Mead  algorithm  is  terminated  when 
the  maximum  least  squares  function  value  difference  between  any  two  points  in  the  simplex 
is  less  than  10-8.  The  algorithm  converges  in  283  iterations  (for  details  see  [BB04b]). 

The  parameters  e^o,  <r  and  the  two  pressure  coefficients  k0 0  of  and  kt  of  r  are 
difficult  for  the  inverse  problem  to  identify.  On  the  other  hand  we  observed  that  the 
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qo 

Iter 

o 

ks 

J5 

l|VJs||L2 

FC 

0.95q* 

33 

80.0999 

48.0607 

4.7019xl0”9 

1.4560xl0-8 

0.0204 

65 

1.05q* 

32 

80.1 

48.06 

9.5079xl0-9 

7.1205x  10-27 

0.0096 

66 

0.5q* 

59 

80.1001 

48.0603 

8.4746xl0-9 

5.5171xl0-9 

0.0130 

112 

Table  2:  Parameter  estimation  of  es,o  and  ks  for  different  initial  simplices  with  no  added 
noise. 


parameters  eS;o  and  its  pressure  coefficient  ns  as  well  as  the  parameter  tq  are  relatively  well 
estimated.  Since  the  inverse  problem  permits  ready  identification  of  only  three  parameters 
in  the  absence  of  noise,  we  do  not  expect  to  see  any  improvement  when  noise  is  added  to 
the  data!  Thus  we  will  henceforth  concentrate  on  the  identification  of  eS;0)  ro  and  ks. 

We  will  first  attempt  the  identification  of  eSio  and  ks.  In  a  first  test,  we  fix  all  the  other 
5  parameters  (including  To)  at  the  true  values  and  attempt  to  identify  eS)o  and  ks.  Table 
2  presents  the  results  for  this  test  using  several  different  initial  simplices.  As  can  be  seen 
from  Table  2,  with  the  other  5  parameters  fixed  at  their  true  values,  we  do  not  have  any 
difficulty  in  identifying  eSjo  and  ns  accurately.  In  Table  2,  FC  denotes  the  function  count, 
i.e. ,  the  number  of  times  the  least  squares  objective  function  is  evaluated. 

The  algorithm  converges  when  the  difference  between  function  values  is  less  than  10 ~8. 
For  the  second  case  where  the  initial  simplex  has  points  with  5%-10%  larger  values  than 
the  true  values,  we  obtain  an  unusually  small  least  squares  function  value  of  the  order  of 
HT27.  We  believe  this  to  be  an  accidental  artifact  of  the  algorithm. 

In  Figures  7  and  8  we  plot  the  least  squares  objective  function  for  values  of  eSio  and  ks 
that  differ  by  0%-10%  from  the  true  values.  The  other  five  parameters  are  kept  fixed  at 
their  true  values.  From  Figures  7  and  8,  we  determine  that  our  inverse  problem  is  more 
sensitive  to  eSio  than  to  the  pressure  coefficient  ns.  This  is  expected  as  the  static  relative 
permittivity,  eS;o  is  more  influential  in  the  system  dynamics  than  the  pressure  coefficients, 
and  hence  is  also  more  important  in  identifying  and  characterizing  the  material. 

As  explained  in  Section  5,  in  a  more  realistic  situation  we  will  not  have  access  to  the 
true  values  of  any  of  the  parameters.  Hence,  we  now  fix  five  of  the  seven  parameters  at 
values  that  have  a  relative  error  from  the  true  values  of  about  10%-50%.  We  use  the  values 

6.0, 

10.0  x  10-12, 

1.5  x  10"5,  (79) 

4.8, 

5.0  x  10~13, 

in  our  simulations.  Such  choices  for  the  parameters  are  typical  in  testing  algorithms  [BK89]. 
We  use  different  initial  guesses  for  eSio  and  ks.  As  discussed  earlier,  noise  is  introduced  into 
the  data.  We  first  add  relative  random  noise  to  our  data,  and  we  use  the  fixed  values  (79) 
in  the  inverse  problem.  The  results  are  tabulated  in  Table  3. 

From  Table  3  we  see  that  using  the  inverse  problem  we  are  able  to  estimate  eSio  but 
we  are  unable  to  reconstruct  ns  with  good  accuracy.  Table  3  shows  the  final  estimates  of 
e^o  and  ns  for  water.  The  other  five  parameters  are  fixed  at  the  values  given  in  (79)  and 
relative  noise  of  varying  levels  is  added  to  the  data. 


€oo,0  — 

TO  = 

a  = 

ft OO  — 

Kr  = 
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Least  Squares  Objective  Function 
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Figure  7:  The  least  squares  objective  function  for  different  values  of  eS)o  and  ks. 


Least  Squares  Objective  Function  Values 


Figure  8:  Cross-sections  from  Figure  7  of  the  least  squares  objective  function  versus  ks 
(left)  and  versus  eSjo  (right).  Note  the  difference  in  scales  between  the  two  figures. 
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To  understand  why  ks  is  not  recovered,  we  note  that  the  fixed  value  of  tq  in  (79)  has 
a  relative  error  of  23%  from  its  true  value.  Since  the  reflections  from  the  acoustic  region 
depend  strongly  on  the  value  of  To,  a  relative  error  of  23%  is  simply  too  high  to  get  a  good 
estimate  of  the  pressure  coefficients.  To  test  this  we  lower  the  value  of  to  to  8.91  x  1CP12. 
So  instead  of  the  values  in  (79)  we  will  use  the  values 

6.0, 

8.91  x  10”12, 

1.5  x  10"5,  (80) 

4.8, 

5.0  x  10~13. 

This  value  of  to  has  a  relative  error  of  10%  from  its  true  value  Tq  .  With  this  new  fixed  value 
for  to  we  attempt  the  inverse  problem  again  to  reconstruct  eSio  and  ks.  These  results  are 
presented  in  Table  4. 

From  Table  4  we  observe  that  the  estimated  value  of  ks  is  closer  to  the  true  value  than  in 
the  previous  test.  Thus  we  can  obtain  accurate  estimates  of  ks  by  choosing  values  of  to  that 
are  closer  to  its  true  value.  We  can  conclude  that  if  the  value  of  tq  is  not  accurately  known, 
then  the  identification  of  ks  is  difficult.  The  identification  of  e^o  on  the  other  hand  does  not 
appear  to  be  sensitive  to  the  fixed  value  of  to  that  is  chosen.  The  combined  identification 
of  both  e<j;o  and  ks  appears  to  be  quite  insensitive  to  the  level  of  relative  noise  that  is  added 
to  the  data. 

We  next  attempt  the  inverse  problem  with  simulated  data  to  which  we  add  constant 
variance  noise.  We  would  like  to  compare  the  results  of  this  case  with  the  results  obtained 
from  the  inverse  problem  with  relative  random  noise  in  the  simulated  data.  We  use  the 
notion  of  signal  to  noise  ratio  to  make  such  comparisons. 


€oo,0  — 

TO  = 

a  = 

^OO  — 

Kr  = 


6.2  Signal  to  noise  ratio 

In  analog  and  digital  communications,  the  signal  to  noise  ratio,  denoted  as  SNR,  is  a  measure 
of  the  signal  strength  (usually  given  in  decibels  dB)  relative  to  the  background  noise.  There 
are  many  equivalent  ways  of  defining  the  SNR  [JN84] .  Since  we  assume  the  mean  of  our 
noise  to  be  zero,  we  define  the  SNR  as 

m=1]  fc' -  dB,  (81) 

Ef=1var  (Ok-E*k)J 


%  RN 

Iter 

ks 

J5 

l|VJ*||L2 

FC 

0.0 

37 

80.3636 

66.0486 

8.5719x  10_y 

2.5594 

0.0076 

73 

1.0 

41 

80.3638 

66.0422 

3.2939x  10”9 

25.2195 

0.0098 

80 

3.0 

40 

80.3640 

66.0280 

8.6023x  10~9 

208.1146 

0.0080 

77 

5.0 

40 

80.3643 

66.0157 

6.3563x  10-9 

574.4429 

0.0060 

79 

Table  3:  Parameter  estimation  of  eSjo  and  ns  for  0%,  1%,  3%  and  5%  relative  noise. 
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%  RN 

Iter 

Cs,  0 

ks 

|J£+i  -  Jf?l 

l|VJ5'||L2 

FC 

0.0 

38 

80.2165 

55.2945 

1.5487xl0-9 

0.4369 

0.0013 

73 

1.0 

39 

80.2169 

55.2914 

4.5901xl0~9 

23.2556 

0.0075 

76 

5.0 

40 

80.2178 

55.2729 

8.2103xl0~9 

573.1135 

0.0079 

78 

10.0 

39 

80.2191 

55.2509 

2.5361xl0"9 

2292.24730 

0.00701 

76 

Table  4:  Parameter  estimation  of  eSio  and  ks  for  0%  -  10%  added  relative 
random  noise.  The  values  of  the  other  five  parameters  are  fixed  at  values 
given  in  (80). 


where  dB  denotes  decibels.  If  the  signal  strength  is  equal  to  the  variance  of  the  background 
noise,  i.e., 

M  M 

^|^.|2  =  ^var(Ofc-^.),  (82) 

fc=l  k= 1 

then  SNR  =  0,  and  the  signal  is  almost  unreadable,  as  the  noise  level  severely  competes 
with  it.  Ideally,  the  strength  of  the  signal  is  greater  than  the  variance  of  the  noise,  so  that 
the  SNR  is  positive.  The  SNR  ratio  can  sometimes  be  increased  by  providing  the  source 
with  a  higher  level  of  signal  output  power  if  necessary.  In  wireless  systems,  it  is  always 
important  to  optimize  the  performance  of  the  transmitting  and  receiving  antennas. 

Using  the  notation  presented  in  Section  5,  the  SNR  in  the  cases  of  relative  random  noise 
SNRr  and  constant  variance  noise  SNRc  are,  respectively,  given  as 


SNRr 


SNRC 


10  log10 


10  log10 


\  Sfc=l  var{l/Ekrlk) 

(  Efeli  \Et\2  ) 

V£f=ivar  (J3m)) 


dB. 


(83) 


To  compare  the  case  of  constant  variance  noise  with  the  case  of  relative  random  noise,  we 
proceed  as  follows. 

1.  As  noted  in  (56)  and  the  discussion  thereafter,  we  set  (3  =  vot.  Setting  the  SNR  in 
the  two  cases  to  be  equal  i.e.,  SNRr  =  SNRc,  implies 


M 

v'Ar(,yEkrik) 

k= i 


Mi/2  a2 


a  = 


(84) 


2.  Thus  for  a  given  value  of  i/,  we  have  (3 
as 

SNRr  =  SNRC  =  10  log10 


av  and  the  corresponding  SNR  is  calculated 
£M)=101ogw(U).  (85) 
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SNR  (dB) 

V 

%  RN 

(3  =  ua 

66 

0.0005 

0.1 

0.0118 

46 

0.005 

1.0 

0.1175 

32 

0.025 

5.0 

0.5877 

26 

0.05 

10.0 

1.1754 

Table  5:  Comparison  of  relative  random  noise  and  constant  variance  noise 


Hence,  for  each  value  of  v  there  is  a  corresponding  value  of  /3  for  which  the  SNR  for  relative 
random  noise  is  equivalent  to  the  SNR  for  constant  variance  noise.  In  Table  5  we  present 
the  values  of  the  SNR,  f3,  and  u  that  correspond  to  each  other  in  the  manner  described 
above.  For  this  test  problem  a  =  23.5077. 

We  next  attempt  to  identify  the  two  parameters  es,o  and  ks  by  adding  constant  variance 
noise  as  given  in  (56)  to  our  data  and  using  the  fixed  values  (80).  The  results  for  the 
corresponding  inverse  problem  are  presented  in  Table  6.  The  results  in  this  case  are  quite 
different  from  the  case  of  relative  noise  as  shown  in  Table  4.  We  can  see  that  the  estimation 
of  ks  becomes  worse  as  the  level  of  noise  increases.  The  estimation  of  eSjo  °n  the  other  hand 
is  relatively  stable,  though  it  is  not  as  accurate  as  the  estimate  obtained  in  the  case  where 
relative  noise  is  added  to  the  data.  The  estimates  of  eSio  also  become  worse  as  the  level  of 
constant  variance  noise  increases.  In  Table  6,  v  corresponds  to  the  scaling  factor  in  (55). 
The  constant  variance  case  and  the  relative  noise  case  for  the  same  value  of  u  have  the  same 
signal  to  noise  ratio. 

We  now  calculate  confidence  intervals  for  the  estimates  that  we  obtained  in  Table  6. 
The  procedure  for  calculating  these  intervals  was  outlined  in  Section  5.1.  These  intervals 
are  presented  in  Table  7  for  varying  levels  of  noise.  From  our  earlier  discussion  about 
confidence  intervals,  we  note  that  the  size  of  the  intervals  is  in  direct  proportion  to  the 
level  of  uncertainty  that  we  have  about  the  estimates  that  we  have  obtained.  As  we  can  see 
from  Table  7  we  have  tighter  intervals  for  the  estimates  of  eS)o  than  we  do  for  the  estimates 
of  ks,  which  is  expected  as  our  problem  is  more  sensitive  to  eSio  than  it  is  to  the  pressure 
coefficient  ks.  We  also  notice  that  the  intervals  become  larger  as  the  level  of  noise  that  is 
added  to  the  data  is  increased.  These  intervals  correspond  to  a  95%  confidence  level. 


%  RN 

Iter 

<bs,o 

ks 

Mn+l  —  Jo  1 

I|VJ5'||L2 

FC 

0.0 

38 

80.2165 

55.2945 

1.5487xl0-9 

0.4369 

0.0013 

73 

1.0 

41 

80.1654 

56.4531 

3.1244  x  10-9 

27.0778 

0.007561 

82 

5.0 

39 

79.9621 

61.1779 

7.8005  x  10~9 

667.7927 

0.0111 

79 

10.0 

44 

79.7084 

67.3235 

3.5288  x  10_1° 

2670.5115 

0.007417 

85 

Table  6:  Parameter  estimation  of  eSjo  and  ks  for  0%-10%  constant  variance 
noise.  The  other  five  parameters  are  fixed  at  the  values  given  in  (80) 
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RN  =  0.0% 

£S,0 

ks 

(8.02165  ±  1.8804  x  10“a)  x  10 
(5.52945  ±  1.6448  x  10~2)  x  10 

RN  =  1.0% 

£s,  0 

ks 

(8.01654  ±  1.4577  x  10~2)  x  10 
(5.64531  ±  1.2738  x  10"1)  x  10 

RN  =  5.0% 

£s,  0 

ks 

(7.99621  ±  6.8206  x  10“2)  x  10 
(6.11779  ±  5.9447  x  10"1)  x  10 

RN  =  10.0% 

0 

ks 

(7.97084  ±  1.2723  x  10-1)  x  10 
(6.73235  ±  1.1076)  x  10 

Table  7:  Confidence  intervals  for  the  parameter 
estimation  of  es,o  and  ks  for  constant  variance 
noise  with  the  Nelder-Mead  algorithm. 


Plots  of  the  least  squares  objective  function  vs.  (eS;o>fts)>  with  the  other  parameters 
fixed  at  either  the  true  values  or  those  in  (80),  and  1%  RN  (constant  variance)  in  the  data, 
reveal  that  indeed  the  minimum  value  for  ns  is  no  longer  near  the  true  value  (see  [BB04b] 
for  details). 

In  Figure  9  we  plot  the  least  squares  residual  versus  time  in  the  case  when  10%  relative 
random  noise  is  added  to  the  data.  A  magnified  plot  of  the  least  squares  residual  versus  time 
in  the  interval  where  the  acoustic  reflection  is  observed  is  plotted  in  Figure  10.  Finally  we 
plot  the  least  squares  residual  versus  the  absolute  value  of  the  electric  field  data  in  Figure 
11  again  in  the  case  where  10%  relative  random  noise  is  added  to  the  data.  We  notice  in 
this  figure  that  the  amount  of  noise  increases  as  the  absolute  value  of  the  electric  field  data 
magnitude  increases.  This  is  in  direct  contrast  to  the  case  when  constant  variance  noise  is 
added  to  the  simulated  data.  We  plot  analogous  plots  of  the  least  squares  residual  for  10% 
added  constant  variance  noise  in  Figures  12  and  13.  Comparing  Figures  12  and  9  we  note 
the  lack  of  any  patterns  in  the  residual  in  Figure  12  for  the  case  of  added  constant  variance 
noise,  whereas  in  Figure  9,  the  residual  follows  the  peaks  in  the  simulated  data  used  in  the 
inverse  problem.  Again  in  Figure  11  we  note  the  fan  like  structure  of  the  residual  with 
the  residual  increasing  with  the  absolute  value  of  the  electric  field.  In  the  case  of  added 
constant  variance  noise,  however,  we  again  note  the  lack  of  a  pattern  in  the  residual  as  seen 
in  Figure  13.  We  also  observe  that  there  are  many  points  of  the  residual  on  the  line  Ex  =  0. 
This  is  because  in  the  simulated  data,  the  electric  field  magnitude  is  almost  zero  most  of 
the  time.  We  compare  the  simulated  data  with  added  relative  random  noise  and  added 
constant  variance  noise  in  Figures  14-17.  We  notice  in  Figure  17  that  the  acoustic  reflection 
is  drowned  by  the  constant  variance  noise  with  v  =  0.05  (10%  RN),  which  is  not  the  case 
with  the  relative  random  noise  data  for  the  same  value  of  u.  Thus  the  identification  of  ks 
is  more  difficult  with  constant  variance  noise  added  to  the  simulated  data. 

We  next  attempt  the  identification  of  three  parameters,  namely  e^Oj  to  and  ks.  The 
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Figure  9:  Plot  of  the  least  squares  residual  versus  time  in  the  case  of  10%  relative  random 
noise  added  to  the  data. 


Figure  10:  A  magnified  plot  of  the  least  squares  residual  versus  time  in  the  interval  where 
the  acoustic  reflection  is  observed,  for  the  case  of  10%  relative  random  noise  added  to  the 
data. 
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Figure  11:  Plot  of  the  least  squares  residual  versus  the  absolute  value  of  the  electric  field 
data  for  the  case  of  10%  relative  random  noise  added  to  the  data. 
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Figure  12:  Plot  of  the  least  squares  residual  versus  time  for  the  case  of  constant  variance 
noise  added  to  the  simulated  data. 
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Figure  13:  Plot  of  the  least  squares  residual  versus  the  absolute  value  of  the  electric  field 
for  the  case  of  constant  variance  noise  added  to  the  simulated  data. 


other  four  parameters,  600,0,  ^00  and  kt  are  fixed  at  values  given  in  (80).  Tables  8  and 

9  present  the  final  estimates  and  the  details  of  the  corresponding  Nelder-Mead  simulation. 
We  observe  that  the  final  estimates  for  both  to  and  ks  start  to  increase  away  from  their  true 
values  as  the  level  of  constant  variance  noise  is  increased.  The  plots  of  the  least  squares 
objective  function  presented  before  provide  an  explanation  for  this. 

In  the  case  of  the  estimation  of  three  parameters  we  observe  that  the  estimates  of  ks  are 
more  accurate  when  tq  is  allowed  to  vary  as  opposed  to  the  case  of  two  parameter  estimation 
when  To  is  fixed  at  a  particular  value.  Varying  to,  however,  does  not  seem  to  make  much 
difference  in  estimating  eSjo-  Also,  comparing  the  intervals  in  Tables  7  and  10  we  find  that 
we  obtain  tighter  confidence  intervals  for  the  three  parameter  case  as  opposed  to  the  two 


%  RN 

Iter 

es,  0 

f0(xlO"12) 

ks 

0.0 

78 

80.1069 

8.1543 

48.2026 

1.0 

90 

80.0745 

8.2242 

49.8265 

5.0 

99 

79.9315 

8.5015 

56.7926 

10.0 

83 

79.7087 

8.8299 

66.2803 

Table  8:  Parameter  estimation  of  eS)o  and  To  and  ks  for 
constant  variance  noise,  corresponding  to  varying  lev¬ 
els  of  relative  noise  with  the  Nelder-Mead  algorithm. 
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Figure  14:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added.  The  value  of  u  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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Figure  15:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  input  signal  at  the  center  of  the  antenna.  The  value 
of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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Figure  16:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  first  reflection  from  the  air-Debye  interface.  The 
value  of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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Figure  17:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  reflection  from  the  acoustic  wave  region.  The  value 
of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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%  RN 

l|VJ5'||L2 

FC 

|J£+i-J£l 

0.0 

0.008942 

140 

8.02867xl0-9 

4.4391  xl0~5 

1.0 

0.001414 

161 

8.5522x  10~9 

26.7071 

5.0 

0.007697 

181 

6.7874x  10~9 

667.6441 

10.0 

0.001432 

154 

6.2983x  10~9 

2670.5056 

Table  9:  Parameter  estimation  of  eS)o  and  to  and  ks  for 
constant  variance  noise,  corresponding  to  varying  levels  of 
relative  noise  with  the  Nelder-Mead  algorithm. 


RN 

0.0% 

es,  o 

TO 

ks 

(8.01069  ±  1.8738  x  10~5)  x  10 
(8.1543  ±  2.3240  x  10"4)  x  10”12 
(4.82026  ±  2.6832  x  10“4)  x  10 

RN 

1.0% 

<±,0 

TO 

ks 

(8.00745  ±  1.4327  x  10“2)  x  10 
(8.2242  ±  1.7739  x  10"1)  x  10'12 
(4.98265  ±  2.0555  x  10“4)  x  10 

RN 

5.0% 

es,  o 

TO 

ks 

(7.99315  ±  6.7761  x  10~2)  x  10 
(8.5015  ±  8.3067  x  10”1)  x  10“12 
(5.67926  ±  9.7751  x  10~4)  x  10 

RN 

10.0% 

^,0 

TO 

ks 

(7.97087  ±  1.2738  x  10"1)  x  10 
(8.8299  ±  1.5286)  x  10"12 
(6.62803  ±  1.8366)  x  10 

Table  10:  Confidence  intervals  for  the  parameter  es¬ 
timation  of  es,ChTo  and  ks  for  constant  variance  noise 
with  the  Nelder-Mead  algorithm. 
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parameter  case,  for  the  same  level  of  added  noise.  The  reason  for  this  again  is  most  likely 
the  dependence  of  the  acoustic  reflections  on  the  value  of  To-  Thus,  the  algorithm  for  the 
inverse  problem  can  choose  tq  appropriately  so  that  the  value  of  ks  is  minimized.  As  before, 
estimates  for  all  three  parameters  become  worse  as  the  level  of  noise  added  is  increased  and 
the  corresponding  confidence  intervals  become  larger  as  seen  in  Table  10. 

6.3  Simulation  Results:  The  Levenberg-Marquardt  method 

We  now  repeat  our  inverse  problem  calculations  with  a  different  least  squares  optimization 
technique.  The  Nelder  Mead  method  presented  in  the  previous  section  is  a  gradient  free 
method  based  on  the  use  of  simplices.  We  will  now  use  a  gradient  based  method,  namely 
the  Levenberg-Marquardt  method  to  solve  our  inverse  problem.  For  the  overdeternrined 
least  squares  objective  function 

1  M  1 

J*(q)  =  -j^2\Ex{ tk,xc,Zc,<i)  -Osk |2  =  -R(q)TR(q),  (86) 

z  fc=i 

Levenberg-Marquardt,  which  is  a  trust  region  method,  adds  a  regularization  parameter 
vlm  >  0  to  the  approximate  Hessian  of  the  Gauss-Newton  method  to  obtain  a  new  estimate 
Qc  =  Qc  +  s,  where 

s  =  -  ( vlmI  +  i?,(qc)Ti?'(qc))'1  R'(qc)TR(qc),  (87) 

with  qc,  the  current  estimate  and  I  the  l  x  l  identity  matrix,  where  l  is  the  number  of 
parameters  that  we  are  attempting  to  estimate.  The  matrix  (vlmI  +  7?/(qc)r7?/(qc))  is 
positive  definite.  The  parameter  vlm  is  called  the  Levenberg-Marquardt  parameter.  For 
details  of  this  implementation  of  the  Levenberg-Marquardt  method  we  refer  the  reader  to 
[Kel99] . 

As  with  Nelder-Mead,  we  first  attempt  the  identification  of  all  the  seven  parameters, 
then  we  identify  e^o  and  ns  and  finally  we  attempt  the  identification  of  eSio>  ro  and  ns- 

Table  11  presents  the  results  for  the  inverse  problem  of  identifying  the  seven  parameters. 
We  observe  how  close  the  final  estimates  of  the  parameters  here  are  with  those  computed  by 
Nelder-Mead.  We  can  again  make  conclusions  similar  to  those  for  the  Nelder-Mead  results. 
Only  eS)o>  to  and  ns  seem  to  be  identifiable  with  some  accuracy,  whereas  the  problem  appears 
to  be  insensitive  to  changes  in  the  other  four  parameters. 

We  next  attempt  the  identification  of  eS)o  and  ks  with  simulated  data  in  which  constant 
variance  noise  is  added.  The  values  of  the  remaining  five  parameters  are  fixed  at  those 
in  (80).  Table  12  presents  the  final  estimates  and  the  details  of  the  corresponding  LM 
simulation.  The  confidence  intervals  for  these  estimates  are  presented  in  Table  13.  Again, 
we  can  make  remarks  and  observations  similar  to  those  for  results  of  the  Nelder-Mead 
method. 

Finally  we  attempt  to  identify  the  three  parameters  e^O;  ro  and  ks.  Table  14  presents 
the  final  estimates  and  details  of  the  algorithm.  Again,  as  was  observed  in  the  results  of 
the  Nelder-Mead  runs,  the  estimates  for  ks  are  better  when  the  value  of  To  is  kept  variable 
as  opposed  to  fixing  the  value  of  to- 

In  Table  15  we  present  the  confidence  intervals  for  the  three  parameter  estimation  prob¬ 
lem. 
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Parameter 

True  Values 

Final  Estimates 

80.1 

80.0917 

^oo,0 

5.5 

5.2039 

TO 

8.1e  —  12 

8.0666  x  10-12 

a 

1.0  x  10~5 

9.5  x  10~6 

ks 

48.06 

48.0028 

Hoo 

4.4 

4.1786 

Kt 

4.05  x  10“13 

3.8462  x  10"13 

Table  11:  Final  estimates  for  all  seven  parame¬ 
ters  for  water  using  Levenberg-Marquardt.  The 
initial  parameter  set  has  values  that  are  5% 
smaller  than  the  true  values. 


%  RN 

Iter 

es,o 

Pvs 

J5 

l|VJ*||L2 

0.0 

18 

80.2166 

55.2951 

0.4369 

0.0006526 

1.0 

18 

80.1654 

56.4524 

27.0778 

0.00076 

5.0 

18 

79.9621 

61.1765 

667.7927 

0.001385 

10.0 

19 

79.7084 

67.3217 

2670.5115 

0.001417 

Table  12:  Parameter  estimation  of  eSjo  and  ks  for  con¬ 
stant  variance  noise,  corresponding  to  varying  levels 
of  relative  noise  with  the  Levenberg-Marquardt  algo¬ 
rithm. 


RN  =  0.0% 

ks 

(8.02166  ±  1.8804  x  10“3)  x  10 
(5.52951  ±  1.6448  x  10“2)  x  10 

RN  =  1.0% 

£s,  0 

ks 

(8.01654  ±  1.4577  x  10“2)  x  10 
(5.64524  ±  1.2738  x  10"1)  x  10 

RN  =  5.0% 

£s,  0 

ks 

(7.99621  ±  6.8206  x  10”2)  x  10 
(6.11765  ±  5.9447  x  10"1)  x  10 

RN  =  10.0% 

£s,0 

ks 

(7.97084  ±  1.2723  x  10_1)  x  10 
(6.73217  ±  1.1076)  x  10 

Table  13:  Confidence  intervals  for  the  parameter  esti¬ 
mation  of  eSio  and  ks  for  constant  variance  noise  with 
the  Levenberg-Marquardt  algorithm. 


40 


%  RN 

Iter 

Cs,  0 

r0(xlO-i2) 

ks 

Ja 

l|VJ*||L2 

0.0 

17 

80.1083 

8.1562 

48.2113 

4.4383xl0-5 

0.0009637 

1.0 

17 

80.0765 

8.2247 

49.8060 

26.7071 

0.001371 

5.0 

19 

79.9320 

8.4983 

56.6654 

667.6441 

0.0008044 

10.0 

20 

79.7085 

8.8330 

66.3476 

2670.5056 

0.001268 

Table  14:  Parameter  estimation  of  eSjo  and  To  and  ks  for  constant  vari¬ 
ance  noise,  corresponding  to  varying  levels  of  relative  noise  with  the 
Levenberg-Marquardt  algorithm. 


RN 

0.0% 

es,  o 

A) 

ks 

(8.01083  ±  1.8741  x  10“5)  x  10 
(8.1562  ±  2.3234  x  10"4)  x  10~12 
(4.82113  ±  2.6830  x  10~4)  x  10 

RN 

1.0% 

£s,  o 
TO 

ks 

(8.00765  ±  1.4334  x  10~2)  x  10 
(8.2247  ±  1.7742  x  lO”1)  x  10“12 
(4.98060  ±  2.0560  x  10"1)  x  10 

RN 

5.0% 

<^£1,0 

TO 

ks 

(7.99320  ±  6.7841  x  10”2)  x  10 
(8.4983  ±  8.3159  x  10"1)  x  10”12 
(5.66654  ±  9.7843  x  10_1)  x  10 

RN 

10.0% 

£s,  o 
TO 

ks 

(7.97085  ±  1.2736  x  10'1)  x  10 
(8.8330  ±  1.5277)  x  10"12 
(6.63476  ±  1.8360)  x  10 

Table  15:  Confidence  intervals  for  the  parameter  estima¬ 
tion  of  e^oTo  and  ks  f°r  constant  variance  noise  with  the 
Levenberg-Marquardt  algorithm. 


7  Sensitivity  of  the  parameter  estimation  problem  to  the  re¬ 
laxation  constant  and  the  interrogating  frequency 

In  [BB04b]  we  have  also  considered  a  second  parameter  estimation  test  problem  with  a 
different  interrogating  frequency  for  a  Debye  medium  with  a  relaxation  time  larger  than 
that  considered  in  the  first  test  problem.  This  was  done  in  order  to  study  the  effect  of  the 
relaxation  time  r  and  the  interrogation  frequency  uic  on  the  ability  of  electromagnetic  pulse 
to  penetrate  the  dielectric  material.  As  shown  in  [BBLOO]  a  very  small  relaxation  time,  such 
as  that  for  water  Tq  =  8.1  x  10~12,  will  make  the  material  appear  comparatively  “hard”  in 
the  sense  that  it  will  allow  much  less  of  the  signal  to  penetrate  the  air-Debye  interface  at 
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z  =  z\  for  a  sufficiently  high  interrogation  frequency  like  coc  =  ir  x  1010  rad/sec.  This  small 
transmitted  wave  can  only  generate  very  little,  if  any,  reflection  at  the  second  interface  as 
it  enters  the  region  containing  the  pressure  wave.  This  demonstrates  the  importance  of  the 
choice  of  the  carrier  frequency  ujc  for  effective  interrogation.  In  [BBLOO]  the  authors  have 
demonstrated  the  strong  influence  of  the  dielectric  relaxation  time  on  an  appropriate  choice 
of  the  interrogating  frequency  for  successfully  penetrating  the  material.  In  our  second  test 
example  (reported  on  in  detail  in  [BB04b])  we  considered  a  parameter  estimation  problem 
for  a  Debye  medium  with  the  relaxation  time  Tg  =  3.16  x  10~8  (as  used  in  [BBLOO,  ABR02]), 
and  an  interrogation  frequency  ojc  =  ir  x  1010  rad/sec. 

As  observed  in  Section  4.1,  the  product  loctq  strongly  determines  the  parameters  in  the 
problem  that  can  be  identified  accurately.  We  observe  that  for  the  choices  of  values  of  tq 
and  loc  in  this  second  test  problem,  the  parameters  that  dominate  in  the  term  e*  in  (51) 
are  no  longer  the  same  as  those  of  the  first  test  problem,  i.e. ,  eS)o  and  ks .  The  parameters 
whose  magnitude  dominates  the  quantity  e*  are  e^o  and  the  pressure  coefficient  k0 0.  This 
observation,  as  was  done  for  the  first  test  problem,  can  also  be  deduced  from  an  analysis  of 
(31).  In  the  simulations  for  this  test  problem  the  outgoing  and  reflected  radiation  will  be 
dominated  by  frequencies  near  the  center  frequency  5.0  GHz.  Thus,  e*  will  be  dominated 
by  frequencies  near  5.0  GHz.  In  this  problem  lotq  «  0(1O2),  and  with  eo  =  8  x  10“12  we 
have  (compare  with  the  estimates  of  (52)) 
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Thus  we  can  see  that  e*  will  be  most  sensitive  to  the  infinite  frequency  permittivity  e^o 
and  the  effects  of  eSio  and  To  and  a  will  not  be  as  pronounced.  Also,  this  implies  that  e* 
will  be  more  sensitive  to  the  pressure  coefficient  Kqo  than  to  the  coefficients  ns  and  kt. 

As  shown  in  [BB04b],  this  second  parameter  estimation  problem  is  very  different  than 
the  one  that  we  have  considered  in  this  paper.  We  summarize  the  essential  differences 
between  these  two  problems  below. 


1.  The  parameters  that  can  be  accurately  estimated  in  the  two  test  cases  are  different.  In 
the  first  example  of  the  Debye  medium  that  is  similar  to  water  we  found  that  eS)o  and 
the  corresponding  pressure  coefficient  ns  can  be  accurately  identified.  In  the  second 
example,  however,  it  is  e^o  and  the  corresponding  pressure  coefficient  k0 0  that  can 
be  accurately  identified.  The  parameters  that  can  be  identified  in  an  inverse  problem 
clearly  depend  on  the  magnitude  of  the  product  loctq- 

2.  In  the  case  of  the  second  test  problem,  which  has  a  larger  relaxation  time,  we  see 
presence  of  local  minima  in  the  objective  function.  The  distance  between  the  local 
minima  increases  as  the  interrogating  frequency  is  decreased,  since  then  the  product 
ujct o  decreases.  Thus,  by  choosing  a  lower  interrogating  frequency  the  local  minima 
could  be  pushed  out  of  the  domain  of  interest.  In  contrast  we  do  not  see  the  presence 
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of  local  minima  in  the  first  parameter  estimation  problem  which  has  a  much  smaller 
value  of  to- 

3.  The  behavior  of  the  Nelder-Mead  algorithm  for  the  parameter  estimation  problem 
for  the  first  Debye  medium  is  markedly  different  than  for  the  second  example.  In 
the  example  with  local  minima  present,  the  Nelder-Mead  algorithm  converges  to  a 
local  minimum  instead  of  the  global  minimum.  We  do  not  see  this  behavior  with  the 
Levenberg-Marquardt  algorithm  for  the  same  example.  In  contrast,  in  the  parame¬ 
ter  estimation  problem  for  water  both  Nelder-Mead  as  well  as  Levenberg-Marquardt 
converge  to  the  global  minimum  and  demonstrate  similar  behavior. 

8  Conclusions  and  future  directions 

In  this  paper  we  have  presented  an  electromagnetic  interrogation  technique  for  identifying 
dielectric  parameters  of  a  Debye  medium  by  using  acoustic  pressure  waves  as  virtual  reflec¬ 
tors  for  the  incident  electromagnetic  pulse.  We  considered  the  two  dimensional  TE  mode 
of  Maxwell’s  equations  which  incorporates  the  pressure  dependence  in  the  Debye  medium 
via  a  model  for  orientational  polarization.  As  a  first  approximation  we  assumed  that  the 
dielectric  parameters  eSio,  £oo,0  and  the  relaxation  tq  are  affine  functions  of  pressure.  We 
have  used  a  finite  difference  time  domain  scheme  to  discretize  our  electromagnetic/acoustic 
model  and  to  compute  simulated  data. 

We  formulated  a  least  squares  based  inverse  problem  for  the  identification  of  the  material 
parameters  as  well  as  the  coefficients  of  pressure  in  the  Debye  model.  We  used  two  different 
methods  to  solve  the  inverse  problem;  namely  Levenberg-Marquardt  and  Nelder-Mead.  The 
inverse  problem  calculations  indicate  that  the  coefficients  that  can  be  identified  depend  on 
the  order  of  magnitude  of  ujtq.  Thus  in  the  first  test  problem  we  were  able  to  identify 
eSio  and  the  corresponding  pressure  coefficient  ns,  and  in  the  second  test  problem  we  were 
able  to  identify  e^o  and  the  corresponding  pressure  coefficient  It  is  not  possible  to 
accurately  identify  all  the  seven  parameters  that  describe  the  Debye  medium  coupled  with 
the  model  for  pressure  dependence  in  either  of  the  examples  that  we  have  considered  in  this 
paper. 

We  have  also  computed  confidence  intervals  for  all  the  estimates  obtained  using  statis¬ 
tical  error  analysis  based  on  the  assumption  of  constant  variance  noise.  We  demonstrated 
that  the  intervals  become  larger  as  the  level  of  added  noise  is  increased. 

There  are  still  many  questions  to  be  investigated.  We  have  used  a  simple  linear  model 
for  the  pressure  term  in  Maxwell’s  equations.  Perhaps  more  importantly,  we  have  assumed 
that  the  effect  of  the  electromagnetic  field  on  the  acoustic  wave  is  negligible.  A  dynamic 
coupled  model  of  the  electromagnetic/acoustic  interaction  is  the  topic  of  future  efforts. 
We  will  also  consider  models  in  which  the  dielectric  parameters  are  nonlinear  functions 
of  pressure  in  attempts  to  understand  better  the  dynamical  pressure  dependence  of  the 
dielectric  parameters. 
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